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ABSTRACT 

Supersonic  flow  past  oscillating  flat  plate  cascades 
with  supersonic  leading-edge  locus  is  analysed  using  a  lin- 
earized method  of  characteristics  valid  for  arbitrary  fre- 
quencies and  an  elementary  analytical  theory  valid  only  for 
low  frequencies  of  oscillation.   These  two  methods  are  ex- 
tensions of  previous  work  by  Teipel  and  Sauer  for  the  single 
airfoil  in  an  unbounded  supersonic  flow  to  the  case  of  air- 
foils oscillating  in  cascade.   Included  is  the  determination 
of  pressure  distribution  and  both  a  two-degree-of -freedom 
(bending  and  torsi'on)  flutter  analysis  and  a  single-degree- 
of-freedom  (torsion)  flutter  analysis.   Numerically  deter- 
mined flutter  boundaries  are  presented  for  various  primary 
parameters  such  as,  Mach  number,  solidity,  stagger  angle, 
density  ratio,  structural  damping  coefficient,  and  elastic 
axis  position.   Also,  results  are  presented  for  the  related 
problem  of  supersonic  wind  tunnel  interference  (including 
the  effect  of  tunnel  porosity)  and  airfoil-airfoil  interfer- 
ence . 
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I.   INTRODUCTION 

In  order  to  design  an  effective  supersonic  compressor,  an 
analysis  of  the  complex  unsteady  flow  phenomena  in  this  speed 
regime  is  necessary.   Since  the  infinite  cascade  has  long  been 
used  by  the  engineer  to  model  two-dimensional  compressor 
flows,  the  study  of  oscillating  supersonic  cascades  is  import- 
ant in  predicting  the  flutter  characteristics  or  dynamic  res- 
ponse of  these  compressors.   Moreover,  the  basic  understanding 
of  the  flow  gained  by  such  a  study  can  lead  to  better  design 
criteria . 

Two  basic  cascade  configurations  can  be  distinguished  at 
the  outset;  i.e.,  cascades  with  either  subsonic  or  supersonic 
leading-edge  locus.   Although  the  former  is  the  case  of  pri- 
mary interest  in  current  research,  in  this  paper  the  case  of 
supersonic  leading-edge  locus  is  treated  in  order  to  form  a 
basis  for  extension  to  the  more  complicated  case  of  subsonic 
leading-edge  locus.   This  method  of  approach  also  has  the 
advantage  that  the  problem  of  oscillatory  wind  tunnel  inter- 
ference, already  analysed  in  previous  work,  is  contained  in 
the  theory  as  a  special  case. 

The  interference  problem  of  linearized  supersonic  flow 
past  airfoils  oscillating  between  solid  wind  tunnel  walls 
was  considered  by  Miles  (1956)  who  derived  a  solution  using 
Laplace  transform  techniques.   Drake  (1956)  treated  this  case 
for  wind  tunnels  with  free  jet  boundaries  and  later  Drake 
(1957)  gave  a  solution  to  this  interference  problem  for 
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porous-wall  wind  tunnels  also  using  Laplace  transform  methods. 
In  an  extension  of  Miles'  work,  Lane  C1957)  presented  a  solu- 
tion for  supersonic  flow  past  oscillating  cascades  with  super- 
sonic leading-edge  locus  again  using  Laplace  transform  tech- 
niques.  Further  work  in  this  area  was  done  by  Hamamoto  (1962). 
Using  Teipel's  (1962)  method  of  characteristics  approach, 
Platzer  and  Pierce  (1970)  made  an  analysis  of  oscillatory 
supersonic  wind  tunnel  interference  and  with  the  help  of  the  . 
high  speed  computer  were  able  to  predict  the  pressure  distri- 
bution along  an  airfoil  oscillating  with  arbitrary  frequency 
in  solid,  free  jet,  or  porous-wall  wind  tunnels.   Platzer 
(1971)  will  present  an  elementary  analysis  of  porous-wall 
wind  tunnel  interference  effects  which  generalizes  Sauer's 
(1950)  solution  for  a  single  airfoil  oscillating  at  low  fre- 
quency in  an  unbounded  linearized  supersonic  flow.   Platzer 
and  Chalkley  (1972)  further  extended  this  solution  to  form  an 
elementary  theory  that  together  with  a  method  of  characteris- 
tics procedure  (based  on  Platzer  and  Pierce,  1970)  is  used  to 
analyse  the  supersonic  flow  past  oscillating  flat  plate  cas- 
cades having  supersonic  axial  velocities  but  otherwise  arbi- 
trary stagger  angle. 

In  this  thesis  the  detailed  description  of  both  the  ele- 
mentary theory  and  the  method  of  characteristics  that  was  not 
possible  in  Platzer  and  Chalkley  (1972)  is  given.   With  the 
aid  of  the  high  speed  computer,  the  method  of  characteristics 
procedure  is  used  to  predict  the  pressure  distributions  and 
aerodynamic  forces  and  pitching  moments  on  a  typical  cascade 


blade.   Both  a  two-degree-of-f reedom  (torsion  and  bending) 
flutter  analysis  and  a  one-degree-of -freedom  (torsion)  flut- 
ter analysis  are  presented,  with  numerical  results  given  for 
the  case  of  the  torsional  flutter.   The  elementary  theory 
(valid  only  for  low  frequencies  of  oscillation)  is  used  to 
provide  a  convenient  check  of  the  characteristics  approach 
as  well  as  to  easily  predict  a  number  of  dynamic  instability 
boundaries  for  a  slowly  oscillating  airfoil  subjected  to 
interference  from  supersonic  wind  tunnel  walls  in  one  case 
and  from  a  larger  airfoil  in  close  proximity  in  yet  another 
case . 
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II.   PROBLEM  FORMULATION 

The  cascade  is  considered  to  be  one  of  two-dimensional 
flat  plates  with  each  blade  performing  equal  low  amplitude, 
simple  harmonic  oscillations  of  the  same  mode  with  the  same 
interblade  phase  angle.   As  shown  in  Figure  1,  it  is  aligned 
in  the  x-y  coordinate  system  such  that  one  blade  lies  along 
the  x-axis  with  its  leading-edge  at  the  origin. 
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Figure  1 

The  cascade  is  further  considered  to  have  supersonic 
leading-edge  locus  and,  except  for  this  restriction,  arbi- 
trary stagger,  such  that, 

tan  g<  cot  «  (1) 

With  this  assumption,  the  flow  between  any  two  adjacent 
blades  can  be  used  to  describe  the  flow  through  the  entire 
cascade.   Since  disturbances  cannot  travel  upstream  of  the 
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Mach  lines  from  the  leading-edge  of  a  blade,  interference 
from  other  than  adjacent  blades  is  not  possible  in  the  super- 
sonic flow  between  them.   Hence,  only  the  flow  between  the 
blade  along  the  x-axis  and  the  one  above  it  is  considered. 
This  flow  is  assumed  to  be  the  non-viscous  flow  of  a  perfect 
gas  governed  by  the  continuity  equation, 

££  +  P  c|i  +  fe  =  o  (2) 

Dt  .   r   dx    dy 
the  Euler  equations, 

dF  +  p  ai    =  °  (3a) 


5X  +  1  |£   .  0 
Dt    p  dy 


(3b) 


and  the  energy  equation, 


§|-0  (4) 

Because  of  the  assumption  of  small  amplitude  oscillations 
of  the  cascade  blades,  all  flow  quantities  are  considered  to 
be  small  perturbations  linearly  superimposed  on  freestream 
quantities.   The  velocities  are  written  as 

u  =  u   +  u  (5a) 

CO  ' 

c  =  c   +  c  (5b) 

CO  N      ' 

while  the  pressure  and  density  perturbations  are  written  as 

Ap  «  p  -  p^  (6a) 

Ap  -  p  -  p^  (6b) 

With  the  assumption  of  isentropic  flow,  the  local  velo- 
city of  sound  can  be  defined  as, 

dp 
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further, 


=  constant 


(8) 


,Y 


Taking  the  total  differential, 

— Z   dP  "  YP   v+i   dp  =  0 
P1  P1 


C9) 


or 


hence , 


dp       p 


c"2  =  y   £ 

P 


(10) 


(11) 


The  flow  boundary  condition  along  the  blade  surfaces  [as 
shown  in  Ch.  5,  Bisplinghof f ,  Ashley,  Halfman  (1955)  for  a 
simple  airfoil]  is  as  follows:   If  the  surface  over  which  the 
fluid  flows  has  the  equation, 

F  (x,  y,  t)  =  0  (12) 

particles  in  contact  with  the  surface  must  have  the  same  nor- 
mal velocity  as  the  surface.   Stated  differently,  the  rate  of 
change  of  F  is  zero  when  the  motion  of  a  particular  fluid 
element  is  followed  along  the  surface,  or 

DF 


Dt 


=  0 


(13) 


If  an  arbitrary,  thin,  two-dimensional  cascade  blade  is 
considered,  the  equation  of  the  upper  surface  can  be  written, 
Fu(x,y,t)  =  y  -  yu  (x,t)  =  0  (14) 

where  y  =  0  is  the  mean  camber  line  and  y   is  the  distance 

*  *  u 

from  the  mean  camber  line  to  the  upper  surface.   Likewise, 
the  equation  of  the  lower  surface  can  be  written, 
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FL  (x,y,t)  =  y  -  yL  Cx,t)  -  0 


(15) 


At  y  -  y 


u 


DF 


u 


3y      3y 


and  at  y  =  y 


Dt 

L 

DF. 
Dt 


at 


3t 


-  u 


-  u 


3x 


Hi 
3x 


+  v  -  0 


(16) 


+  v  =  0 


(17) 


since 


3y     3yT 

J  u     J  L 


=  1 


3y     3y 
Thus,  the  normal  flow  velocity  can  be  written, 


(18) 


3y 


v  = 


v  = 


u 


3y 


3t 

9y, 


+  u 


+  u 


u 


3x 

3y, 


at 


at 


y  -  y„(x,t) 


u 


y  =  yL(x,t) 


(19a) 


(19b) 


3t       3x 

By  applying  the  assumption  of  linear  perturbations  on  the 
freestream,  the  normal  velocity  equations  become 


3y 


v  = 


V  = 


u 


3y 


3t 

Hi 

3t 


+  u 


u 


00  3x 

ay, 


+  u 


3x 


+  u 


+  u 


ay, 

3x 

Hi 

3x 


at  y  -  yu(x,t) 


at  y  =  yL(x,t) 


(20a) 


(20b) 


Since  a  thin  blade  is  being  considered  y   and  yT  are  small, 

u       L 

thus  the  equations  become,  cancelling  the  higher  order  terms, 


3y 


v  = 


v  - 


u 


3y 


3t 

3y, 


+  u 


u 


•  3x 

3y. 


at  y  =  y  (x,t) 


u 


+  u 


at  y  -  yL(x,t) 


(21a) 
(21b) 


3t   '  "»  dt 

This  normal  velocity  can  be  further  written  as  a  Taylor 
series  expansion  of  the  normal  flow  velocity  at  y  =  0.   Thus, 
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v(x,yu,t)  -  v(x,0+,t)  +  yu  8y<*>°  'C) 


.2 

2! 


+  Ik    a2v(x,o,,t)  + 
ay2 


(22a) 


and 


v(x,yL,t)  =  v(x,0  ,t)  +  yL 


3vU,0",t) 


3y 


y£   32v(x,0~,t) 
+  TT     3y2        + 


(22b) 


Again  with  the  assumption  of  a  thin  blade  and  small  linear 
perturbations  on  the  freestream:  y  ,  y  ,  and  v  are  small  and 
higher  order  terms  may  again  be  cancelled  from  the  equations 
leaving : 


v(x,y  ,t)  =  v(x,0  ,t) 


u 


v(x,yL,t)  =  v(x,0",t) 


Thus, 


3y 


u 


3y 


v(x,y,t)  =  y£—  +  u 
v(x,y,t)  =  j^-   +  u 


u 


3x 

3y, 


at  y  =  0 


+ 


°°  3x 


at  y  =  0* 


(23a) 
(23b) 

(24a) 
(24b) 


With  the  assumption  of  simple  harmonic  motion  of  the  blade, 


yu(x,t)  =  yu(x)  e 


yL(x,t)   =  yL(x)  e 


iut 


and  the  normal  velocity  equations  can  then  be  written, 

3y 
v(x,y,t)  =  [iwy  (x)  +  uffl  j^- J  e      at  y  =  0 

3y 
v(x,y,t)  =  [iwyL(x)  +  u^  j^-]    e5^*   at  y  =  0" 


(25a) 
(25b) 

(26a) 

(26b) 


Since  the  flow  is  identical  between  any  two  adjacent 
blades  , 
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3y 

J  e        at  y  =  d     (27) 


v(x,y,t)  =  [iwyL(x)  +  u^ 


In  plunge,  the  lower  blade  oscillation  about  y  ■  0  is, 

iu)t 


hence, 


v(x,y,t)  =  [-iwho]  e 


y  =  -  h   e' 

4  o 


Iw  t 


(28) 


at  y  =  0     (29) 


while  the  upper  blade  oscillation  about  y  =  d  is, 

.    i((ot+5) 
y  =  -  h   e 
o 

where  6  is  the  interblade  phase  angle.   Hence, 

v(x,y,t)  =  [a)  sin  6  h   -  ioo  cos  6  h  ]  e 

at  y  =  d  . 


(30) 


(31) 


In  pitch,  the  lower  blade  oscillation  about  y  =  0  is, 

1(0 1 


y  =  -  9  (x-x  )  e 
o     o 


hence , 

v(x,y,t)  =  [-u 


0   -  iu)6   (x  -  x  )  ]  e 

oo    o  O  O 


i(ot 


(32) 


(33) 


at  y  =  0  ,  while  the  upper  blade  oscillation  about  y  =  d  is, 

i(u)t+<5) 


y  =  -0   (x  -  B  -  x  )  e 
'      o  o 


(34) 


hence , 


v(x,y,t)  =  [-0   {u   cos  6  -U)  sin  6  (x-B-x  )  } 


10)t 


-  16   {u   sin  6  +  0)  cos  6  (x-B-x  )}]  e    (35) 
o    °°  o 

at  y  =  d  . 

The  pressure  distribution  on  the  blades  can  be  written  in 
terms  of  the  local  sonic  velocity  as  follows:   In  terms  of 
perturbation  quantities,  Eq.  (11)  can  be  written, 
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1     Pco  +  AP 

Poo  M  x   A£> 


c2  +  2c  C   =  Y  -    r—  (36) 


cancelling  the  higher  order  c2  term.   Since  — "  is  less  than 

°  °  p 

r  00 

one,  Eq.  (36)  can  be  expanded  in  the  geometric  series, 

1  =    !    .   M+   A£^    -1-4*  (37) 


(1+    A£)  Poo  Poo  Poo 

Poo 

neglecting  the  higher  order  terms.   Substituting  this  in 

Eq.  (36)  gives, 

1        c"   A 

2c  c   =  Y  -±-    (1 -  ££■)  Ap  (38) 

ooip^      y   Ap 


Since  Ap  and  Ap  are  small, 


Ap     d_p_      -2 
dp 


(39) 


Substituting  this  in  Eq.  (38)  gives  after  some  algebraic 
manipulation, 

2   c  c   ( |L )  =  -4-  Ap  (40) 


Y-l     -  \  +   2c         Poo 

(y-1) c 

»  I     *   oo 

Since  2c/(y-l)c   is  small,  x may  be  expanded  in 

i    +  — ±£ 

(Y-Dc 

x   I  '       oo 

the    geometric    series, 

1  -    1     -    TT^V-      +     lT-4f-]2    -     ....     -    1    -  2C 


-,  2c  (Y-l)c  '(Y-l)c    '  (Y-l)c 

-f-     — x  '  '      oo  v,/oo  v  '  '      oo 

(Y-D^ 

(41) 

neglecting  higher  order  terms.   Hence,  Eq .  (40)  can  be  writ- 
ten, neglecting  higher  order  terms, 

2 
P  "  P~  =  7"  P   c   (c  -  c  )  (42) 

■     ■  oo    Y  — 1       °°         °° 
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Finally,  any  non-dimensional  quantities  so  stated  in  the 
paper  are  made  such  with  reference  to  the  blade  chord  length 
and  the  freestream  velocity  unless  otherwise  defined. 


18 


III.   KETHOD  OF  CHARACTERISTICS 

In  the  method  of  characteristics  it  is  desirable  to  as- 
certain the  coordinate  system  across  which  all  possible  dis- 
continuities in  flow  properties  may  occur.   Along  this 
coordinate  system  the  equations  of  motion  of  the  flow  field 
can  then  be  treated  as  ordinary  differential  equations  that 
are  solvable  by  classical  or  numerical  techniques  (e.g., 
finite  differences).   To  obtain  the  equations  of  this  coor- 
dinate system  (the  characteristic  directions)  in  the  (x,y) 
plane,  the  equations  of  motion  are  written  in  terms  of  the 
arbitrary  intersecting  coordinates, 

C  =  5  (x,y)  (43a) 

and 

n  =  n  (x,y)  (43b) 

If  the  first  derivatives  of  the  dependent  variables,  u,  v, 
and  p  (or  c)  with  respect  to  E,    are  made  indeterminate  across 
lines  of  ri  =  constant,  and  the  first  derivatives  of  the  de- 
pendent variables  with  respect  to  r)  are  made  indeterminate 
across  lines  of  E,    =    constant,  then  any  possible  discontinui- 
ties in  the  first  derivatives  will  occur  across  these  lines. 
These  lines  are  then  the  characteristics  and  their  equations 
are  obtained  in  the  evaluation  of  the  inde terminacies . 

Consider,  first,  two-dimensional,  steady  flow.   The  gov- 
erning equations  of  motion  are  the  continuity  equation, 

M£El  +  i|£vi.  =  o  (44) 

ox        9y 
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the  Euler  equations, 

-  3u      3u    1  3p  // c  n 
u  -5—  +v-r—  +  —  t^k  0                            (45a) 

3x      3y    p  3x 

-  3v  ,    3v    1  9p    ~  //c:u\ 
u  ■=—  +  v  -s—  +  —  -r^-  =  0  (45b) 

3x      3y    p  3y 

and  the  energy  equation, 

1 

u|^+v^=  0  (46) 

3x     3y 

Along  a  streamline  this  equation  may  be  written, 

-  3p      3p_„r-3p      3p,  //-7N 

U  Vs1  +  V  -r*-  =  c2[u  -r^-  +  V  -r21]  (4  7) 

3x      3y   c    3x      3y 

Substituting  in  Eq.  (44),  the  continuity  equation  becomes, 

r  3u    3v,      1   ,-  3p      3p  n  /  /  on 

P  [■5—  +  T~]  =  -  -   [u  ^  +  v  -r-^-]  (48) 

r   3x    3y       -2     3x      dy 

In  terms  of  the  arbitrary  coordinates,  Eqs.  (43),  the 

continuity  equation  becomes, 

p5      ur    +    p£      v>.  +  -       [u£      +   vC    ]    Pr 
x      £  y%*2  x  y         £ 

=    -pn      u      -    pn      v      -—      [uri      +vn]p  (49) 

K  x     n  y     n       -2  x  y     rn 

In  like  manner,  the  Euler  equations  become, 

p[u£      +    v£    ]    u      +    £    p^.    =    -p[uri      +   vti    ]    u    -    p    n  (50a) 

p[u£      +   vf    ]    v      +    E    P^    =    -p[uf|      +   vn    ]    v    -    p    n  (50b) 

^x    ^y   £    y  5       x    y   n   n  y 

Eq.  (49)  and  Eqs.  (50)  form  a  system  of  three  simultaneous 
equations  in  u  ,  v  ,  and  p  .   Solving  for  p   by  Cramer's 

5-55  5 

rule  gives  a  ratio  of  two  determinants.   The  denominator  is 
the  determinant  of  the  matrix  of  coefficients  of  ur,  Vr,  and 
Pr  formed  by  Eq.  (49)  and  Eqs.  (50),  while  the  numerator  is 
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this  determinant  with  the  right  hand  side  of  the  equations 
substituted  for  the  coefficients  of  p  . 

Since  across  lines  of  n  =  constant,  p   is  indeterminate, 
both  the  numerator  and  the  denominator  must  equal  zero.   Set- 
ting the  denominator  equal  to  zero  and  expanding  the  deter- 
minant gives, 


[u?x  +  vM  [(u2-  I2)    g*  +  2uv  55   +(v2-  c"2)^]  =  0 


"x^y 


(51) 


The  solutions  to  this  equation  give  the  equations  of  all 
three  characteristics  in  the  physical  (x,y)  plane.   These  are, 


di  =  v  m    t 
dx    u 


4*  =  t. 


an  £ 


(52a) 


.an  (C  ±  a)  (52b) 

where  £  is  the  angle  the  streamline  makes  with  the  x-axis  and 
a  is  the  Mach  angle.   All  three  characteristic  directions  are 
obtained  from  Eq .  (51)  because  F,    and  r)  are  arbitrary  and  in- 
terchangeable:  If  Eq.  (49)  and  Eqs.  (50)  were  written  such 
that  derivatives  with  respect  to  r\    are  the  unknowns,  making 
p   indeterminate  across  lines  of  F    =  constant  would  result  in 

n 

solutions  identical  to  Eqs.  (52),  a  priori.   Eq.  (52a)  is 
merely  the  equation  of  a  streamline  while  Eqs.  (52b)  are  the 
equations  of  the  right  and  left-running  Mach  lines.   Tradi- 
tionally, Eqs.  (52b)  are  the  equations  of  F,    =  constant  and 
H  =  constant. 

The  compatibility  relation  which  relates  changes  in  u,  v, 
and  p  along  the  (£,n)  characteristics  can  be  obtained  by 
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setting  the  numerator  of  p   equal  to  zero.   Expanding  the 

determinant  and  solving  the  resultant  equation  yields  after 

some  algebraic  manipulation, 

3u      3v  r  1        3p    n  /r,x 

v  -t—  -  u  -r—  +  —  cot  a  Tr*-   =  0  (53) 

8n     3n   p      3n 

This  derivation  is  shown  in  detail  in  Chapter  8  of  Oswatitsch 
(1956).   These  relations  must  be  satisfied  along  the  (£,r|) 
characteristics . 

The  physical  significance  of  these  characteristics  is 
that  they  act  as  boundaries  across  which  any  discontinuities 
will  occur.   In  this  flow  field  where  streamlines  and  Mach 
lines  are  the  characteristics,  a  discontinuity  is  readily 
apparent  across  the  Mach  line  emanating  from  the  leading-edge 
of  a  slender  airfoil  or  across  a  slip  plane  (streamline) 
emanating  from  a  Mach  stem.   Characteristics  can  further  be 
described  as  information  carriers  in  that  along  them  all  dis- 
turbances propagate  and  all  flow  quantities  can  be  determined 
given  their  values  in  the  domain  of  dependence.   It  should 
also  be  noted  that  the  characteristic  directions  are  com- 
pletely independent  of  the  type  of  coordinate  system  used  to 
describe  the  flow.   They  are  determined  solely  by  the  physical 
nature  of  the  flow  field  as  expressed  in  the  equations  of 
motion. 

The  problem  of  two-dimensional,  unsteady  flow  over  a  flat 
plate  as  shown  by  Teipel  (1962)  can  be  considered  in  the  same 
manner  as  before.   The  governing  equations  of  motion  [Eq.  (2), 
Eq.  (3),  and  Eq.  (4)]  are  written  in  an  arbitrary  coordinate 
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system  across  which  variation  of  flow  properties  are  indeter- 
minate.  The  evaluation  of  these  indeterminacies  will  again 
give  the  characteristic  directions  of  these  coordinates. 

The- arbitrary  coordinates  are  given  in  Eqs.  (48).  Eq.(2), 
Eq.  (3),  and  Eq.  (4)  can  be  rewritten  in  a  form  more  suitable 
for  transformation  by  using  the  local  sonic  velocity  as  an 
unknown.   Assuming  that  air  is  a  perfect  gas,  Eq.  (11)  can 
be  written, 

(54) 


c2  -  Y  —  T 


where  <%_     is  the  universal  gas  constant  and  ^  is  the  molecular 

weight.   Taking  the  total  differential  and  dividing  by  c 

gives, 

a  ~  /Oat 

(55) 


2  il        A  il 


Substituting  Eq.  (54)  into  this  gives, 


2  l£  -  g*  (56) 

c 

In  like  manner,  the  total  differential  of  Eq.  (11)  is, 


2c  dc   =  y 


£  .  i2.  _  v  £  .  d£ 


p   p 


p   p 


and  from  Eq .  (11)  again, 


dc    dp    dp 


(57) 


(58) 


From  the  First  Law  of  Thermodynamics  and  the  definition  of 
entropy, 


Tds   =  du  -  *-   di 


(59) 
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Since  du  -  c   •  dt,  this  may  be  written  as, 


Tds  «  c  dT  -  RT 

v         p 


da 


(60) 
C61) 


where  R  is  the  gas  constant  (c  -c  )  and  c   and  c   are  the 

p   v        p       v 

specific  heats  of  the  gas  at  constant  pressure  and  volume, 
respectively.  Dividing  Eq.  (61)  by  Tc  ,  substituting  from 
Eq.  (56)  and  Eq.  (59),  and  then  taking  the  substantial  deri- 


vative of 


-S  with  respect  to  time  gives, 


I     £i=2Ti     2£  -    (y-1)    1     52.  =    o 
c       Dt  Y    -       Dt  VY       }    p       Dt 


(62) 


since  entropy  is  held  constant.   This  can  be  written  substi- 
tuting from  Eq.  (11), 


Dp     2    -  Dc 

— *-  =  pc  — 


(63) 


Dt    y_1  *""  Dt 

Introducing  the  small  perturbation  assumption  and  substitut- 
ing from  Eq.  (7)  and  Eq.  (63),  Eq.  (2)  can  be  written,  can- 
celling higher  order  terms, 

8u  ,  8v 


2   8c     2      8c 

=-  tt—  +  — r  u   t: —  +  C  -z —  +  c  -X —  =  0 

y-1  8t    y-1   °°  8x     °°  8x     °°  8y 


(64) 


This  is  the  form  of  the  continuity  equation  that  will  be  used 
in  the  coordinate  transformation.   Again  dividing  Eq .  (61)  by 
Tc  ,  substituting  from  Eq.  (56)  and  Eq .  (58),  and  then  taking 
the  partial  derivative  with  respect  to  x  (holding  entropy 
constant)  gives, 


1  ip_  _ 


-  8c 


(65) 


p  8x    Y""1     ^x 
Substituting  in  the  first  Euler  equation,  Eq.  (3)  becomes, 
cancelling  higher  order  terms, 
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3t     «  3x   Y-l   «  3x 
In  like  manner,  the  second  Euler  equation,  Eq.  (4),  becomes, 
3v       3v  ,   2      3c"    .  ,A_* 

"5T  +  u   -r—  +  r-  C   77-  =  0  (67) 

3t     °°  3x   y-1   »  3y 
The  continuity  equation  and  the  Euler  equations  thus  form  a 
system  of  three  simultaneous  equations  in  u,  v,  and  c.   With 
assumption  of  simple  harmonic  motion,  Teipel  (1962)  intro- 
duced his  amplitude  functions, 

U(x,y)eiU)t 


u    -    u 

oo 

u 

oo 

1 

V 

/m2-i 

u 

00 

2         1 

c    - 

c 

V(x,y)eit0t  =   _L_    -7  (69) 


C(x,y)eia)t  =^r-^  *—  -- 
Y-1   M2   °°° 

where  U,  V,  and  C  are  complex. 

Substituting  in  Eq.  (64),  Eq.  (66),  and  Eq.  (67)  gives, 

|£  +  /M2-!  |I  +  M2  |£  +  ikM2C  =  0  (71) 

dX     r         dy        dx 

H  +  |£  +  ikU  =  0  (72) 

dx     dx 

|Z  +  ^Lr     |£  +    ikV    =    0  (73) 

9x  /M2-l       9y 

where  k  is  the  reduced  frequency. 

Transforming  this  system  into  the  new  coordinate  system, 
£(x,y)  and  i~[(x,y),  and  considering  now  a  system  of  three 
simultaneous  equations  in  unknowns,  U_,  Vp  and  C_,  the  equa- 
tions may  be  written,  (74) 
Cx  V^   +  cot  a  5yVc+  M2?XC5  =  -Un-  cot  a  Vnny-M2Cnnx  -  ikM2C 

ix  us  +  cx  cc    =  -Vx-Vx  "  iku  (75) 
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^x   V£  +  tan    a   ^y    cc    "   "Vn   ^   -   tm  a.  Cn   ny    -  flOT  ((7160. 

As    before    bV  ,    V^. ,    and    C_    axe   jmaidie    inxle LermXTra=rrFe  aExrxss^ 
lines    of    r\    =    constant.       This    forces    aaxy   d i s c on tLnudcties;  inr 
U,    V,    and    C    to    lie    across    lines    of    ri    =    constant.       Oiree  cojxdi:— 
tion    of    indeterminacy    is    that    the    determinant    of    the^  matrix: 
of    coefficients    of    uV  ,    Vp    and    Cr    in    the   above    egni a^joirsc  be 
zero.       If    this    determinant    is    expavmded   and    then    sae^t:  eaqxradl  to: 
zero,    the    resulting    equation    is, 

-  «*  -  e,  k;  -  «*s*]  =  °  ('7^ 

As  before  the  solutions  to  this  equation  give  the  three;  char- 
acteristic directions  in  the  physical  (x,y)  plane: 

421  =  +  tan  a  (73b.) 

ax    — 

These  coincide  with  the  results  af  the  steady  ELaw.  analy- 
sis when  it  is  remembered  that  the  small  perturbation  assump- 
tion forces  the  streamlines  to  be  parallel  with  the:  x-aocis 

tt  -  o). 

In  order  to  find  the  compatibility  relations  of  the  first 
derivatives  of  U,  V,  and  C  along  the  characteristics,  the 
equations  of  motion  are  written  along  lines  of  E,    =  constant, 
T\   -    constant,  and  along  streamlines,  that  is,  for 
str(x,y)  =  constant, 


|J  -  0  (7?a) 


E    (x,y)  =  constant,      -s±   =    _  (73b-) 
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n(x,y)  =  constant,  -^  =  - 


3x 


/K2-l 


(79c) 


If  an  arbitrary  function  of  x  and  y  in  the  equations  is 

3f 


denoted  by  f  then,  -s-—  holding  str  constant  becomes, 

3f 


v3x'str    3x 


(80) 


3f 


Likewise,  y-  holding  £  constant  becomes, 


(44> 


3f 


+ 


3f 


3x  k     3x     /M2-i    9y 


(81) 


3f 


and,  t: —  holding  n  constant  becomes, 


(44) 


3f 


3f 


3x  n    3x  /M2_1  3y 

Combining  these  equations  and  rewriting, 

3f  ,3f. 

3x  ^3x;str 


3x    2L^3x;?    ^3x;nJ 


(82) 


3f 


3f 


3f. 


3y  =1^2-   i  [(^>C  "  Cft>n] 


(83a) 
(83b) 
(83c) 


Substituting  these  relations  in  the  equations  of  motion  gives, 

1 

2 

3x'£ 


2U3x'£ 


3U. 

3x'n 


t[(-^:),  +  £f)J  +  t(m2-  i)[(  ). 


3V, 
3x'£ 


3V. 
3x'ti 


+  TM2[(^),  +  (S)  J  +  ikM2C  =  0 


3C 

3x-'n 


(84a) 


J    "t£>5  +  Cfx-V  +  |l#5  +  Cff)^  +  ikU  -  0       (84b) 

*(!^  +  €V  +  IICfx^  "  C^V  +  ikV  =  0       (84c) 
Subtracting  Eq.  (84b)  from  Eq.  (84a)  gives  an  equation  that 
is  first  added  to  Eq.  (84c)  and  then  subtracted  from  Eq.  (84c) 
The  result  is  two  compatibility  relations  that  must  be 
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satisfied  along  lines  of  constant  £  and  lines  of  constant  T\. 
These  are, 

cH>r  +  Cl^)'   +  ik  IV  +  -^—  (M2C  -  U)]  =  0       (85a) 

^n  '  ^n  +  ik  IV  -    -1—    (M2C  -  U)J  =  0        (85b) 

ox  n    dx  n         m2-i 

The  third  compatibility  relation  is  obtained  by  considering 
Eq.  (84b)  along  a  streamline, 

<~0  ,.   +  (~)  ,   +  ikU   =0       (86) 
dx  str     dx  str 

If  irrotationality  is  assumed,  the  equation  of  irrotationali- 

ty, 

'SH-I*.  0  (87) 

dy        dx 

may  be  substituted  for  the  second  Euler  equation  in  the  sys- 
tem.  Again  with  the  assumption  of  simple  harmonic  motion,  it 
may  be  written  in  terms  of  Teipel's  amplitude  functions, 

|2-  -  y^-^T   H  .  o  (88) 

8y  3x  v 

The  non-dimensional  system  of  equations  can  then  be  written 

as , 

f£  +  /M2-  1  |l   +  M2  P  +  ikM2C  =  0  (89a) 

dx  9y        3x  v 

|£  +  «  +  iku  .  o  (89b) 

f  "   ^^i=°  (8'c) 

By  transforming  Eqs .  (89)  to  the  arbitrary  coordinate  system, 
£  and  ri ,  as  before,  the  system  can  be  considered  one  of  simul- 
taneous equations  of  U  ,  V  ,  and  C  .   By  making  any  of  these 
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derivatives  indeterminate  the  results  of  setting  the  deter- 
minant of  the  matrix  of  their  coefficients  equal  to  zero  are 
Eqs.  (78),  the  characteristic  directions. 

By  writing  the  equations  of  motion,  as  before,  along 
lines  of  £  =  constant  and  n  -  constant  and  performing  the 
same  algebraic  manipulations,  the  two  compatibility  relations 
can  be  determined: 

(^n +  #r, +  lk  -frx  (u  -  c>  ■  °     (90b> 

The  third  compatibility  relation  is,  as  before, 

(~)  ,   +  (|^)  .   +  ikU  =  0  (86) 

dx  s tr     dx  s tr 

The  equations  of  motion  are  now  reduced  to  a  system  con- 
taining only  derivatives  with  respect  to  x.   It  is,  thus, 
possible  to  solve  these  equations  using  finite  differences 
in  the  following  manner:   Separating  the  real  from  the  imagi- 
nary components,  the  equations  can  be  written, 

'jT'str  +  (3^>str  "  MI  '  «  <91a> 

3U  dC 

(■5-^-)  _   +  (-5-^)  *.   +  klT   =  0  (91b) 

dx   s  tr     9x   s tr      R 

3U        3V  2 

3UR        3VR 


M2-l 
M2 

CD,- 

V 

M2-l 
M2 

M2-l 
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9UX       9V  2 


(UR-  CR)  -  0 


(91f) 


To  write  these  equations  in  finite  difference  form,  a 
computational  molecule  as  shown  in  Figure  2  will  be  used. 

With  the  velocity  values  at  P-,-.,  P-io'  and  P21  known»  aH  tlie 
values  at  P„„  can  be  calculated.   The  distance  Ax  is  arbi- 
trary . 


Figure  2 


If  F  is  an  arbitrary  flow  quantity,  its  partial  deriva- 
tives can  thus  be  written, 

(92a) 

(92b) 

(92c) 

The  values  of  the  flow  quantities  in  the  equations  are  aver- 
ages, such  that, 

(F)str  =  f  (F11  +  F22>  C93a) 


(91) 
v8x'str 

F 
22 

Fll 

Ax 

(21) 
3x  5 

=  - 

F 
J22 

yAx 

F12 

(91) 

= 

F22 

2Ax 
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(F)   -  \   (F12  +  F22)  (93b) 

(F)n  =  |  (F21  t  F22)  (93c) 

With  these  relations  substituted  into  the  equations,  the  sys 
ten  becomes,  after  some  algebraic  manipulation, 

U22R  +  C22R  -  AIU22I  "  K12R  <94a) 

U22I  +  C22I  +  AIU22R  "  K12I  <94b) 

U22R  "  V22R  "  BIU22I  +  BIC22I=  K34R  (94c) 

U22I  -  V22I  +  BIU22R  "  BIC22R=  K34I  <94d) 

U22R  +  V22R  "  BIB22I  +  BIC22I=  K56R  (94e) 

U22I  +  V22I  +  BIU22R  "  BIC22R=  K56I  (94f> 


where  : 


and , 


K12R  -  U11R  +  C11R  +  AIU11I  (95a) 

K12I  "  "ill  +  C11I  -  AIU11R  (95b) 

K34R  "  U12R  "  V12R  +  BI(U12I"  C12I)  <95c) 

K34I  =  U12I  "  V12I  "  VU12R"  C12R)  (95d) 

K56R  =  U21R  +  V21R  +  VU21I-  C21I>  (95e) 

K56I  =  U21I  +  V21I  "  BI(U21R"  C21R>  (95f) 


A  =  J  k  4x  (96a) 

1      M2 
B   »  T-  k  — - —  Ax  (96b) 

A    *         M2-l 


As  shown  in  Teipel  (1962),  solving  these  equations  gives, 
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u 


U-VI).li<K34K+W-BIK12I)+2BI[iCK34I«56I)+BIK12]l] 


22R 


(1-Vi)2  +  C2V2 


C97a) 


U 


a-AIBI)l|cK34l+K56l)+BIK12Rj-2BI[|(K34R+K56R)-BIK12I] 


221 


(1-A^)2  +  (2Br)2 


V2  2R  '"  2  (K5  6R~  K34R) 


V22I    2  (K56l"  K34I) 


C22R  "    K12R   U22R+  AIU22I 


C22I    K12l"  U22I~  AIU22R 


(97b) 
(97c) 
(97d) 
(97e) 
(97f) 


These  are  the  finite  difference  equations  for  a  point  in  the 
general  flow  field. 

At  a  point  on  the  top  of  a  cascade  blade  or  solid  bound- 
ary, the  normal  velocity  V9~  is  prescribed  by  the  movement 
of  the  blade.   This  velocity  was  given  in  the  Problem  Formu- 
lation.  Here  the  computational  molecule  is  shown  in  Figure  3 


Figure  3 
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The  applicable  equations  are, 


U22R  +  C22R  ~  AIU22I    K12R 


U22I  +  C22I  +  A1U22R  "  K12I 


U22R  ~  BIU22I  +  BIC22I  "  K56R~  K34R 


U22I  +  BIU22R  '  BIC22R  ~  K56l"  K34I 


C98a) 
C9  8b) 
(98c) 
(98d) 


where  K- 0  and  Kc,  are  as  before,  but, 
iz        jo 


K     =  V 
*34R     22R 


K     =  V 
^341     221 


(99a) 
(99b) 


as    given   by    the    boundary    conditions.       Solving    for    these 
equations    gives, 

(1-AlV[K56R-K34R-BIK12l]+    2  VK56I-K34I+BIK12r1 


U 


u 


22R 


221' 


(1-AIBI)2  +  (2Bx)2 


(1-AIBI)[K56I-K34I+BIK12R]-2BI[K56R-K34R-BIK12I] 


(1-AjBj)2  +  (2BX)2 


C22R"  K12R  ~  U22R  +  AIU22I 


(100a) 

(100b) 
(100c) 


C22I   K12I    U22I    AIU22R 


(lOOd) 


V22£  and  V_„   are  known  from  the  boundary  conditions.   These 
are  the  finite  difference  equations  for  the  flow  velocities 
at  the  top  of  a  solid  boundary  or  cascade  blade. 

Likewise,  at  the  bottom  of  a  solid  boundary  or  cascade 
blade,  the  normal  velocity,  Y  2o    *<s  prescribed  by  the  movement 
of  the  blade,  and  was  given  in  the  Problem  Formulation.   The 
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computational  molecule,  different  from  that  at  the  top  of  a 
blade,  is  shown  in  Figure  4. 


,y 

pn 

AX 

P22 

P12 

X 

— >• 

Figure  4 


Here,  the  applicable  equations  are, 


U22R  +  C22R  "  AIU22I    K12R 


U22I  +  C22I  "  AIU22R  "  K12I 


U22R  '  BIU22I+  BIC22I"  K56R  +  K34R 

U22I  +  BIU22R"  BIC22R=  K56I  +  K34l 
where  K1 _  and  K  ,  are  as  originally  stated,  but, 


K     =  V 
*56R    V22R 


K     =  V 
*56I    V22I 


(101a) 
(101b) 
(101c) 
(lOld) 


(102a) 
(102b) 


as  given  by  the  boundary  conditions.   Thus,  the  finite  differ- 
ence equations  for  the  flow  velocities  at  the  bottom  of  a 
solid  boundary  or  cascade  blade  are, 

U-AlV  [K56K+K34R-BIK12.I]+2BI  [K56I+K34I+BIK12EJ 


U 


22R 


(1-AIBI)2  +  (21^)* 


(103a) 
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u 


221  C1-AXBX)2  +  C2BX)2 


C22R~  K12R  "  U22R  +  AIU22I 


C22I~  K12I  ~  U22I  "  AIU22R 


(103b) 

Q03c) 
(103d) 


V00   and  V__   are  given  from  the  boundary  condition. 
L  Z.  R         &  Z  i. 

To  obtain  the  conditions  on  the  initial  left-running  Mach 

line,  assume  that  P_  ..  and  P91  as  shown  in  Figure  2  are  just 

in  the  freestream  and  then  let  Ax  shrink  to  zero.   Since 

K_  _  =  Kc/.  =  0  and  A_  and  Bx  go  to  zero  as  Ax  goes  to  zero, 
±z     j  b  1       l 

the  initial  finite  difference  equations  for  £  =  constant  are, 

U22  +  C22  =  °  (104a) 


Thus, 


U22  "  V22  +  V12  -  ° 


u22  +  v22   =  o 


U22    "  V22      C22 


(104b) 
(104c) 

(105) 


Eq.  (90a)  can  then  be  written, 


cg>    -  -  iku 


Integrating  gives, 


U  =  U 


Thus, 


U22R  =  -  V22RCO)COS  ^ 


U22I  =  -V22l(0)cos(k 


exp [-ik 
x=0 


M- 


M2-l 


x] 


W 


(106) 


(107) 


Mz-  1 

2 


x)-  V22I(0)sin  (k  — - 


M' 


x)  (108a) 


>r-  l 


M 


M  -  1 


x)+  V22R.(0)sin(k  — 


M' 


x)   (108b) 


M*-  1 
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V22R  =  -  U22R  (108c) 

V22I  -  -  U22I  (108d) 

C22R  =  -  U22R  C108e) 

C22I  "  "  U22I  <108f> 

For  the  initial  right-running  Mach  line,  P  -  and  P  _  of 

Figure  2  are  just  in  the  freestream  and  Ax  is  brought  to  zero 

Since  K- 0  =  K_,  =  0,  and  Ax  =  B_  =  0, 
12     34  II 

U22  +  C22  =  °  (109a) 

U22  "  V22  =  °  (109b) 


U22  '  U21  +  V22  =  °  (109c) 


Thus, 


U22  =  V22  =  "  C22  (110> 

Eq.  (90b)  can  then  be  written, 


Integrating  gives, 


Thus, 


c|a>   -  -  iku  (in) 


U  =  V(0)  eL~liC  m2_    XJ  (112) 


M2    ,  .„    ,_  .  ,,    M2 


U??t?=  V    (O)cos(k  — - x)+V99T(0)sin(k  — - x)       (113a) 

22R    22R  M2-  1      221  M2-  1 

2  2 

U??t=  V?9T(0)cos(k  — —  x-)-V99_(0)sin(k  — — —  x)      (113b) 
221    221  H2-  1      22R  M2-  1 

V22R  =  U22R  (113c) 

C22R  =-U22R  <113d> 
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V22I  ~  U22X 


^221     U22I 


CH3e) 
(113f) 


The  pressure  distribution  along  the  surfaces  of  the  blades 
can  be  determined  from  Eq.  (42).   Dividing  both  sides  by  p^ 
and  using  Teipel's  amplitude  function  for  C,  the  equation 
can  be  written, 


P(x,y)  =  yM2C(x,y) 


x=0 
x=0' 


(114) 


where 


P(x,y)e 


iut 


P  "  P, 


Thus  the  non-dimensional  pressure  distribution  along  the 
top  surface  and  the  bottom  surface  of  the  cascade  blade  can 
be  computed  in  a  point  by  point  calculation  of  U,  V,  and  C 
along  lines  of  £  =  constant  and  r|  =  constant. 
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IV.   ELEMENTARY  THEORY  FOR  SLOWLY  OSCILLATING  CASCADES 

Sauer's  (1950)  solution  for  an  airfoil  oscillating  at 
low  frequencies  in  an  unbounded  supersonic  flow  can  be  applied 
to  the  flat  plate  cascade  to  form  a  theory  suitable  for  com- 
parison with  the  method  of  characteristics.   As  shown  in  Gar- 
rick  (1957),  the  equations  of  motion  of  supersonic  flow  over 
a  flat  plate  oscillating  at  small  amplitude  can  be  written, 
assuming  irrotational  flow,  as 

cot2a$    -  $  +    2   -  M<$>      -  -   $    =0  (115) 

xx     yy     c    xt    2   tt 

00 

With  the  assumption  of  simple  harmonic  motion, 

$(x,y,t)  =   <Kx,y,k)eikt  (116) 

and  the  potential  equation  becomes, 

cot2a<j>    -  <J>   +  2ikM2d)   -  k2M24>  =  0  (117) 

xx     yy         x 

where   <J>  is  complex. 

If  the  frequency  of  oscillation  is  sufficiently  low,  the 
perturbation  potential  amplitude  may  be  expanded  in  a  Taylor 
series,  such  that 

4>(x,y,k)  =  x(x,y)  +   k   ^    (x,y)  (118) 

neglecting  higher  order  terms.   Eq.  (117)  then  splits  into 
two  simultaneous  equations  since  k  cannot  appear  in  any  rela- 
tion, it  being  independent  of  k.   These  equations  are: 

cot2a  Y    ""  X    ■  0  (119a) 

Axx    ^yy 

cot2a  ip        -  i>        -  -2iM2x  (119b) 

Txx    ryy         Ax 
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Sauer  (1950)  showed  that  a  general  solution  to  these  two 
equations  is, 

X  "  gCZ) 


$   =  h(Z)  -  i 


M 


cos  a 


yg(z) 


Tp  =  h(Z)  +  i 


Z  =  x  -  y  cot  a 
M 


yg(z) 


C120) 
(121) 

(122) 
(123) 


cos  a 

Z  =  x  +  y  cot  a  (124) 

where  g  and  h  are  arbitrary  functions  of  position  and  are 
equal  to  zero  for  z  ^  0  (x,y,  z  and  d  are  non-dimensional). 

To  apply  this  general  solution  to  the  cascade  with  super- 
sonic leading-edge  locus,  the  flow  field  is  divided  into 
separate  zones  as  shown  in  Figure  5 .   The  number  of  these 
zones  is  dependent  on  the  size  of  A.   For  A  >  1,  there  is 
one  zone  along  each  blade;  for  0.5  <  A  ^  1,  there  are  two; 
for  0.33  <  A^  0.5,  there  are  three;  for  0.25  <  A<  0.33, 
there  are  four;  etc. 

If  oscillation  in  pitch  only  is  analysed,  the  boundary 
conditions  on  the  blade  surfaces  can  be  written,  in  non- 
dimensional  form,  from  Eq .  (33)  and  Eq .  (35), 

at  y  =  0         (125) 


<t>y  =  -1  -  ik(x  -  xq) 

<f>   =  ei6[-l  -ik(x  -  x  )] 


at  y  =  d 


(126) 


where  9   is  one. 
o 
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Figure  5 

In  Zone  I  there  is  no  interference  from  the  upper  blade, 
hence,  Sauer's  solution  for  the  single  airfoil  applies.  The 
boundary  condition,  however,  must  be  satisfied  at  y  =  0, 
hence , 

(127a) 
M 


Xy  =  -  cot  a  g0'(Z)  =  -1 


\b      -    -    cot  a  h  '  (Z)  -i 
y  o        cos  a   o 


g^(Z)=-i(x-xo)  (127b) 


This  gives, 


g  (Z)  =  Z  tan  a 
°o 


(128a) 


h  (Z)  =-i  tan  aZ(x  +  =r   tan2a)  (128b) 
o  o   2 


where , 


M 


*  =  g0(Z)+  k{hoCZ)-i  ^f-^y  go(Z)} 


(129) 


In  terms  of  Teipel's  amplitude  functions  the  flow  properties 
can  be  related  to  those  obtained  in  Method  of  Characteristics 

by, 
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U    =       <f>x  (130a) 

V   =    tan    a    4>  (130b) 

C    =    -Hx   +    ik<}>]  (130c) 


In    Zone    I, 


U  =  tan  a  (1  -  ik[x  +  x  tan2a  +  y  cot  a]}  (131a) 

V  =  tan  a  {-1+  ik[x  -  x  +  y(cot  a  +  — - — )]}       (131b) 
y I  o  cos  a   J 

C   =  tan  a  {-1+  ik[x  +  x(tan2a-l)+  2y  cot  a]}       (131c) 

In  Zone  I'  the  perturbation  potential  can  also  be  written 
directly  from  Sauer's  single  airfoil  solution, 

*  ■  8o<*>  +  k  {ho^  +  1  t^Hj  y  s0(?)}  (132) 

where , 

Z=x-A-B+y  cot  a  (133) 

This  ensures  that  Z  <;  0  upstream  of  the  leading-edge  of  the 
upper  blade . 

The  boundary  condition  at  y  =  d  gives, 

Xy  =  cot  a  g^(Z)  =  -  el6  (134a) 

ill   =  cot  a  h'(Z)+  i  — - —  g  (Z)+  i  — —cot  ag' (Z) (134b) 
ry  o        cos  a  &o        cos  a      &o 

hence,  =  -  i(x  -  x   -  B)  e 

o  .  ~ 

gQ  (Z)  -  -  Z  tan  a    e10  (135a) 

h  (Z)  =  i  tan  a  Zei6{x  +  f  tan2a+  — — }       (135b) 
ov  o   2         cos  a        v 

Then, 

UIt=    -tan    a    eXQ{l-ik!x   +(x-B)tan2a   +A-y    cot    a] }  (136a) 

VT,=    -tan    a    eX<5{-lfik[x    -  (x-B)  +  (d-y  )  (cota+   — )]}  (136b) 

i.  o  "  cos    ot 

C    ,=    -tan    a    ei<5  {-1+ik  [x   +(x-B)  (tan2a-l)+2A-2y    cot    a]  }     (136c) 


41 


For  0.5  <  A  <  1.0,  Zones  II  and  II'  will  be  present  in 
addition  to  Zones  I  and  I1.   In  Zone  II  the  perturbation 
potential  is  due  to  the  initial  waves  from  Zones  I  and  I1 
plus  the  reflected  waves  at  Zone  II,  hence, 

<P   -  gQ(z)  +  go(z)  +  g1(z1)+k{ho(z)+ho(z) 

M 


+  h.  (Zn)-i 

1   1     cos  a 


y[g0(z)-  go(Z)-g1(Z1)]}   (137) 


where  Z  and  Z  are  defined  in  Eq.  (122)  and  Eq.  (133),  and, 

(138) 
The  boundary  condition  at  y  =  d  gives, 

16 


Z1    =  x  -  2A  +  y  cot  a 


Xy  =  -  cot  a  g^(Z)+  cot  a  g^(Z)+  cot  a  z[^±)    =-e: 
ipy    =  -  cot  a  h^(Z)+  cot  a  h^(Z)+  cot  a  h-JCZ.) 


M 


1  c^T^  [g0(z)-  80(z)  "  «i<zi>] 


+  i  - 


Md 


cos 


-  cot  a  [gJ(Z)+  g^(Z)+  gjU^] 
±6 


=-i(x-x   -B)e 

hence,  after  some  algebraic  manipulation 

g1(Z1)  =  Z   tan  a 

h  (Z  )  =  -i  tan  a  Z. {x  +  ^  tan2a  +  2    Md  } 
x   x  l   o   z  cos  a 

Then, 

UII=  tan  a{2  -  ei<5  -ik [ (2-ei6) (Xq+  x  tan2a  +  A) 

+  e   (y  cot  a  +  B  tan2a)]} 
VIZ=  tan  a{-e1Q+  ik[el0(x  -  x  +(d-y)(cot  a+  — - — )) 


cos    a' 


-2(d-y) (cot    a   + 


M 


cos    a 


)]} 


CXI=    -tan    ot{2-ei<S  +    iklei<5  (xot(x-B)  ( tan2a-l) 
+    2A    -    2y    cot    a)-    2 (xq+   x(tan2a-l)    +    2A) ] } 


(139a) 


(139b) 

(140a) 
(140b) 


(141a) 


(141b) 


(141c) 
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In  Zone  II1  the  perturbation  potential  is  due  to  the  ini- 
tial waves  from  Zones  I  and  I'  plus  the  reflected  waves  at 
Zones  II  and  II',  hence, 

<J>  -  g0(z)  +  s0(z)+  g1(z1)t  k{hQ(Z)+  ho(z) 

+  VV-  i  Tolb  ^S0(Z)-  *o<Z>+  8l(Zl>])  (1A2) 

where  Z  and  Z  are  defined  in  Eq.  (122)  and  Eq.  (133),  and, 

Z   =  (x  -  A  -  B)-  y  cot  a         (143) 


The  boundary  condition  at  y  =  0  gives, 

X   =  -  cot  a  g^(Z)+  cot  a  g^(Z) 

-  cot  a  g-fC^)  =  -1 

ip      =  -cot  a  h'(Z)+  cot  a  h'(Z)-  cot  a  h'(Z  ) 


-i 


M 


cos  a 


[gQ(Z)-  gQ(Z)+  g1(Z1)]  =  i(x  -  xq) 


(144a) 


(144b) 


hence, 


g-  (Z.  )  =  -  Z-  tan  a  e 


±6 


V  in     \  n  1(5  r         1      2         Md 

h.(Z.)=  i  tan  a  Z.  e   ix   +  -r—  tan  a  +  

11  1       o    2  cos  a 


} 


(145a) 
(145b) 


Then, 

UTT,=  tan  a{l-2e    +  ik[2e   (x  +(x-B)tan2a 
II  o 


+  A)  -  (x  +  x  tan2a  +  y  cot  a)]} 

VTT,=-tan    a{-e1<S+    ik  [  2e1<5  (y  (cot    a   +   — - — )) 
II  w  cos    a 

-    (x    -    x   +   y(cot    a   +   — - — ))]} 
o  ' v  cos    a      J 

CXI,=    -tan    a{    l-2ei<5  +    ik[2ei<5(x   +(x-B) (tan2a-l) 
+    2A)    -    (xq+   x(tan2ct-l)  +    2y    cot    a)]} 


(146a) 


(146b) 


(146c) 
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Because  of  the  algebraic  complexity  of  further  analysis, 
additional  zone  calculations  will  be  made  with  the  assumption 
that  the  blades  are  oscillating  in  phase,  i.e.,  <5  =  0 . 

For  0.33  <  A  ^  0.5  in  addition  to  Zones  I,  I',  II,  and 
II',  Zones  III  and  III'  will  be  present  in  the  flow  field. 
In  Zone  III  the  perturbation  potential  is  due  to  the  initial 
waves  from  Zones  I  and  I'  plus  the  reflected  waves  at  Zones 
II,  II1,  and  III, 
<f>  =  80(Z)+  S0(Z)+  g1(Z1)+  g1(Z1)+  g2(Z2)  + 

k{ho(Z)+  hQ(Z)+  h1(Z1)+  h1(Z1)+  h2(Z2) 


-l 


My 


cos 


-tg0(Z)-  80(Z)"  S^V*  gl(Zl)+  g2(Z2)]}     (147) 


where  Z,  Z,  Z..  ,  Z.     are  defined  as  before,  and, 

Z2  =  x  -  2A  -  y  cot  a  (148) 

The  boundary  condition  at  y  =  0  gives,  after  considerable 
algebraic  manipulation, 


g2  =  Z2  tan  a  (149a) 

Z, 
2    o   2  cos  a 


2      o  MH  , 

h   =  -  tan  a  Z.  {x  +  ^  tan2a  +  2  — -^— }  (149b) 


Then, 

UIII  =  tan  a^1_  ik^x  +(x~2B)tan2a  +  2A+  y  cot  a]}      (150a) 

VTTT  =-tan  a{l-  ik[x  -  x  +  y(cot  a  +  — - — )]}  (150b) 

ill  o       v  cos  a 

C1T1    =-tan  a{l-  ik[x  t (x-2B) (tan2a-l)+4A  +2y  cot  a]}   (150c) 

In  Zone  III'  the  perturbation  potential  is  due  to  initial 
waves  from  Zones  I  and  I'  plus  the  reflected  waves  at  Zones 
II,  II1 ,  and  III1 , 
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<j>    =    gQ(Z)+    gQCZ)+    g1(Z1)+    g^Z^)*    g2(Z2) 
+    k{hQ(Z)+    hQCZ)+    h1(Z1)f    h1CZ1)-f-    h2CZ2) 

-    i    dfV8o<Z>-    *o&-    *1CV  +    glCV~    ^2^2)]}  <151> 

where    Z,    Z,     Z1 ,    and    Z       are    defined    as    before,    and 


Z2=x-B-3A+y    cot    a 


(152) 


The  boundary  condition  at  y  =  d  gives,  after  considerable 

algebraic  manipulation, 

(153a) 
(153b) 


g0(Z0)  =  -  Z9  tan  a 
Z. 


2      2 
h0(Z0)  =  i  tan  a  Z0  {x  +  -r—  tan  a  +  3 


Md 


2  v  2 
Then, 


o   2 


cos  a 


} 


UTxxt=  -  tan  a{l  -  ik[x  +(x-3B) tan2a+  3A  -  y  cot  a]}   (154a) 
III  o 

VTTT,=  -  tan  ct{l-ik[x  - (x-B)  + (d-y ) (co t  a  +  — - — )]}    (154b) 
III'  o  J  cos  a 

Cj   ,=  tan  a{l-ik[x  + (x-3B) ( tan2 a-l)+6A  -  2y  cot  a]}   (154c) 

For  0.25^  A^  0.33,  in  addition  to  Zones  I,  I',  II,  II*, 
III,  and  III1,  Zones  IV  and  IV1  will  be  present  in  the  flow 
field.   In  Zone  IV  the  perturbation  potential  is  due  to  the 
initial  waves  from  Zones  I  and  I1  plus  the  reflected  waves 
at  Zones  II,  II',  III,  III1, 
<t>  =  gQ(Z)+  gQ(Z)+  g1(Z1)+  81(Z1)+  g2(Z2)+  g2(Z2)+  g3(Z3) 

+  k{ho(Z)+  hQ(Z)+  h1(Zi)+  h1(Z1)+  h2(Z2)+  h2(Z2)+  h3(Z3) 
My 


-x 


cos 


-Ig0(z)-go(z)-g1(z1)tg1cz1)+g2(z2)-g2(z2)-g3(z3)j} 

(155) 


where  Z  through  Z_  are  defined  as  before,  and 

Z3  =  x  -  4A  +  y  cot  a 


(156) 
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The  boundary  condition  at  y  =  d  gives  after  considerable 
algebraic  manipulation, 

g3(Z3)  =  Z3  tan  a  (157a) 


Z 


3  ^   2   .  /    Md 


h0(Z.)  =  -i  tan  a  Z.{x  +  -^-   tan2cc  +  4  — -^—}  (157b) 

3   3  3   o   2  cos  a 

Then, 

UIV=   tan  a{l-ik[x  +(x+3B)tan2a  +  3A  +  y  cot  a]}       (158a) 

VT  =  -tan  a{l-ik[x  - (x-B)+ (y-d) (cot  a+  — - — )]}        (158b) 
IV  o         J  cos  a 

C   =  -tan  a{l-ik[x  +(x+3B) ( tan2 a-l)+6A+2y  cot  a]}      (158c) 

In  Zone  IV'  the  perturbation  potential  is  due  to  the  ini' 
tial  waves  from  Zones  I  and  I1  plus  the  reflected  waves  at 
Zones  II,  II1,  III,  III',  and  IV': 

$   =    g0(Z)+    g0(Z)+    S^V*    81(V+    82(Z2)+    82(V  +    83(V 
+k{hQ(Z)+    hQ(Z)+    h1(Z1)+    h1(Z1)+    h2(Z2)+    h2(Z2)+    h3(Z3) 

'-[g^(Z)-g^(Z)-gl  (Zj+g.  (2)+g,(Z0) 


cos    al6ov    '    6ov    '     61V    1'     61V    ]/     &2V    2' 

-g2(Z2)+    g3(Z3)]}  (159) 

where  Z  through  Z_  are  defined  as  before,  and, 

Z=x-B-3A-y  cot  a  (160) 

The  boundary  condition  at  y  =  0  gives, 

g3(z\)  =  -  Z   tan  a  (161a) 


Z 


3  _  2.   .  .    Md 


h,(Z.)  =  i  tan  a  Z_  {x  +  -^-   tan2a  +  3  — l-^—}  (161b) 

3   3  3    o   2  cos  a 

Then,  after  considerable  algebraic  manipulation, 

UT„t=  -  tan  a(l-ik[x  +(x-4B)tan2a  +  4A  -  y  cot  a] >     (162a) 
iv  o 

VTV,=    -    tan    ctU-iklx    -    x    -    y(cot    a   +   — - — )]}  (162b) 

IV  o  v  cos    a 
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C   ,=  tan  a(l-ik[x  - (x-4B) Ctan2 a-l)+  8A  -  2y  cot  a] }   (162c) 

Further  zones  may  be  calculated  for  values  of  A  less  than 
0.25  in  the  same  manner. 

In  the  Method  of  Characteristics  it  was  shown  that  the 
non-dimensional  pressure,  P,  is  a  direct     function  of  C. 
By  integrating  this  pressure  over  the  blade  surface,  the  lift 
and  pitching  moment  on  the  blade  may  be  determined.   The  deri- 
vation of  the  lift  and  moment  equation  is  given  in  the  Flutter 
Analysis.   In  terms  of  Garrick  and  Rubinow's  (1946)  lift  and 
moment  coefficients,  for  the  case  of  in-phase  blade  oscilla- 
tion (6  =  0), 

2    *      +        1+B 
L  =  L   +  i  L,  =  -  =—{f    C(x,0  )dx  -  /     C(x,d  )dx}     (163) 

J       *      k2  o  B 

4  ,  l  + 

M  =  M_  +  i  M.  =  -  — {/  (x-x  )  C(x,0  )dx 

3       4      k2  o 

l+B 
-  /     (x  -  B-  xo)C(x,d  )dx}  (164) 

B 

For  0.25  <  A  <  0.33, 

L3  =  ° 

L4  =  8  ^  tan  a{d2  +  3A2  +  B2(tan2a-1)}  •  (165) 

M0  =  -16  —  tan  a{  A2  -  B2  }  (166) 

3        k2 

M4=  -16  £  tan  a{x  (d2+2A2+B2 tan2a) -8A3 -4d2 A+4AB2 }      (167) 

If  the  cascade  blades  are  allowed  to  oscillate  at  an  ar- 
bitrary interblade  phase  angle,  the  pitching  moment  determi- 
nation can  be  somewhat  simplified  by  integrating  in  the 
following  manner:   For  0.5  ^  A  <  1.0, 


47 


M  =  M-+  i  M,=  -  — {  /     (x-x  )ICT(x,0+)~ 
3. .     *      k2   ° 

-  ei0  CT,(x+B,d  )]dx  t  /   Cx-x  )ICT(x,0T ) 
1  A-B      P    k 

.  «  1  . 

e    CTT(x+B,d~)]dx  +  /  (x-x  )[CTT, (x,0  ) 
11  A+B     °    •LJ" 


-ei<S  CI];(x+B,d")]dx} 


(168) 


This  integration  gives 

M.  =  8  -  tan  a{cos  6 [x  (2-2A)+(A-B) 2 -1-2AB]  -k  B  sin  6 


{2x2  +  x  [2A+2(tan2a-l) ] -4A2 - ( tan2a-l) 
o     o 


(A2+  1  -  ^B2)}  -  x„ 


1 
2 


(169a) 


M.  =  8  f  tan  ct{cos  6[x2(2A-2)+  x  (2A2-  4A  +  2 
A      k  o  o 

+  (A2-  B2-  l)tan2a)-  -j  A3+  2(1-B2)A  -  j 

+  -I  (1-A3)tan2a]  -  2  ^  B  sin  6  [x   -  A] 
3  k  o 


+  x2  -  x   +  ^  +  (^x   -  h    tan2a} 
o     o    3     2  o    3 


(169b) 


In  most  airfoil  theory,  the  non-dimensional  pitching 

moment  is  written  as, 

C   =  9  [c    +  ik  c   ]eikt  (170) 

m     ol  mQ        m'QJ 


Relating  this  to  Garrick  and  Rubinow's  (1946)  moment  coeffi- 


cients gives  (where  6   is  one) , 


1  ,_2 


me       2      3 


mG      2     4 


(171a) 
Cl71b) 
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V.   FLUTTER  ANALYSIS 

Lane  (1956)  showed  that  the  flutter  analysis  of  a  cascade 
can  be  made  by  considering  only  a  single  blade  of  the  cascade 
oscillating  in  the  flow  field.   In  order  to  compute  the  non- 
dimensional  lift  and  moment  acting  on  that  blade,  the  pressure 
difference  over  the  top  and  bottom  surfaces  must  be  integrated 

over  the  blade.   Thus, 

l        + 
L  =  /'  {P(x,0  )  -  P(x,0  )}dx  (172) 

o 

i  + 

M  =  /   (x-  xq){P(x,0  )-  P(x,0  )}dx         (173) 

o 
where  P  is  the  non-dimensional  pressure  as  defined  by  Eq. 
(114). 

Lift  here  is  positive  downward  as  in  Garrick  and  Rubinow 
(1946).   The  non-dimensional  pressure  on  the  lower  surface  of 
this  equivalent  blade  is  that  of  the  lower  surface  of  the 
blade  above  adjusted  for  the  interblade  phase  angle  and 
stagger  such  that, 

P(x,0~)  =  e1<S  P(x+B,d)  (174) 

The  lift  and  moment  on  the  equivalent  blade  can  also  be 

written 

1      +        1+B   -16 
L  =  «f  P(x,0T)dx  -  /     e     P(x,d)dx  (175) 

o  B 

M  =  /  (x-x  )p(x,0T)dx  -  /  (x-B-x  )e  igP(x,d)dx 
o      °  B         ° 

°  (176) 

Separating  these  into  their  real  and  imaginary  parts  and  writ- 
ing the  equations  in  terms  of  the  local  sonic  velocities  gives, 
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,     1  +  2+B 

L  =  yM2{/  C_,(x,0T)dx  -  /[cos  6>CD(x,d)  +  sin<5«CT(x,d)]dx} 


Cl77a) 


L  =  yM2{/  CT(x,0T)dx  -  /Icos  <$«CT(x,d)-sin  6  •  C_  (x,  d)  J  dx} 

B  (177b) 

,  1  +        1+B 

M  =  yH2{f    (x-x  )CD(x,<r)dx  -  /  (x-B-x  ) 

R      o       °   R  b 

•[cos  6  C_(x,d)+  sin  5  •  CT  (x, d) ] dx}  (178a) 

1  .         l+B 

MT=  yM2{/  (x-x  )CT(x,0  )dx  -  /     (x-B-x  ) 
I      ~       o   I  _  o 


•[cos  6-CT(x,d)-  sin  6-C_(x,d) ]dx}  (178b) 

where  C   and  C   are  the  real  and  imaginary  components,  res- 
pectively . 

To  evaluate  these  integrals  numerically,  the  trapezoidal 

rule  is  used : 

1  Ax 

/   f(x)dx  =  ~[f(x)  +  2f(x).  +  2f(x)_+  2f(x)_  +  .. 
2.  o         1  z         J 

o 

+  2f(x)   ?+  2f(x)   .+  f(x)  1        (179) 

n-z        n-l       n 

If  oscillation  only  in  bending  and  torsion  (pitching  and 

plunging)  is  considered  and  Garrick  and  Rubinow's  (1946) 

method  of  expressing  lift  and  pitching  moment  is  used, 

1  2h 

L  -  ~  f  Poo  £  "I   k2e1U3t[— £(LX+  i  L2)+  6o(L3+  i  L4)J    (180) 

c 


2h 

< 


M  =  -  £  9(o    c2u2  k2eia)t[— ^CM,+  i  M0)  +  9  (M,+  i  M ,  )  J    (181) 


toe 
where  k  =  - —  unlike  Garrick-Rubinow ' s ,  h   is  the  vertical 
u  o 

00  w 

equilibrium  position  of  the  elastic  axis,  and  0   is  the 
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equilibrium  angle  of  attack  of  the  blade.   In  dimensional 
form  the  lift  and  pitching  moment  equations  are, 


L  -  c  pm  e 


iwt 


{[Lp+  i  LTJ      +  IL..+  i  LJ       } 
K     *  Plunge   K     x   Pitch 


mode 


mode 


(182) 


~2      .lWt 
oo 


'  =  c'pm  eltut  {[MB+  i  MT] 


+[M  +  i  M  ]      } 
K        Plunge  Pitch 


mode 


mode 


(183) 


Thus, 


L,  +  i  L0  =  - 

1  2        m2,.2 


—  [L   +  i  L  ] 
M'k^Y  h  Plunge 

mode 


(184) 


Separating  the  real  from  the  imaginary  parts  of  the  equation, 


L   =  - 

1        ^t2i.2 


M^k^y  h 


'R   Plunge 
mode 


(185) 


or, 

?       l        .         l  +  B 

Lt= {/  C_(x,0  )dx  -/[cos  6  •CTJ(x,d)+sin  6 

k   h   °  B 


•CT(x,d)]dx}_1         ,  (186a) 

I  Plunge  mode 


Likewise  , 


l+B 


L  =  - 

*■  *  2 


k'  h   o 

o 


{/  C].(x,0+)dx  -  /[cos  6  •  C][(x,d)  -  sin 


•CR  (x,d)]dx}plun 


ge  mode 


(186b) 


l        ,         l  +  B 
L„  =  -2  {/  C_(x,(T)dx  -  /[cos  6  •  CTJ(x,d)+  sin  6 

3      k2  6   o    R  B  R 

o 


•CT(x,d)  Jdx},,.,.  .     , 
I  Pitch  mode 


(187a) 
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L  = 


1.1       +        l+B 

-2  — {/  CT(x,0T)dx  -  /  Icos  <$-CT(x,d) 

k2  9   o    X  B  X 


-  sin  6«C_(x,d)]dx} 


Pitch  mode 
l+B 


M,=  -2 


k'  h   o 
o 


l  ,        l+tf 

{/  (x-x  )C_,(x,0T)dx  -  /   (x-B-x  ) 

O    is. 


B 


•[cos  S«C_(x,d)+  sin  6 • CT (x , d) ] dx} 
K  J. 


Cl87b) 


Plunge  mode  (188a) 


i+B 


M  =  -2 


k^  h   o 


l  .         l^B 

{f    (x-x  )C  (x,0  )dx  -  /  (x-B-x^) 


B 


M  =  -4 


•[cos  6*C  (x,d)-  sin  <$  •  C  (x,d)]dx} 

Plunge  mode  (188b) 

1      1  +        1+B 
{/  (x-x  )C_(x,0  )dx  -  /  (x-B-x  ) 

k2  6   o       °   R  B 

o 


•[cos  6«C  (x,d)+  sin  6«C  (x,d)]dx) 

Pitch  mode   (189a) 
l+B 


,       1  .  l+B 

= {/  (x-x  )C  (x,0  )dx  -  /  (x-B-xq) 


M,=  -4  — 

k2  9   ° 
o 

•[cos  6»C  (x,d)-  sin  6»C  (x,d)]dx> 

Pitch  mode   (189b) 

where  the  integrals  are  computed  using  the  trapezoidal  rule 
and  the  two  different  oscillatory  modes  are  determined  by 
the  boundary  conditions  at  the  blade  surfaces. 

Bisplinghof f ,  Ashley,  and  Halfman  (1955)  define  flutter 
as  "the  dynamic  instability  of  an  elastic  body  in  an  air- 
stream."   In  Chapter  9  they  describe  the  physical  nature  of 
the  flutter  phenomenon.   Basically,  flutter  occurs  when  the 
characteristic  determinant  of  the  equations  of  motion  of  the 
blade  vanishes. 
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If  blade  oscillation  in  two  degrees  of  freedom  (pitch  and 
plunge)  is  considered,  the  equations  of  motion  of  the  cascade 
blade  can  be  written  as, 


Mh  +  SQ9  +  Chh  =  L 


sQh  +  i  e  +  cQe  -  m0 


where  the  total  inertia  force  acting  on  the  blade  is, 

fx  =  -  (Mh  +  s0e) 

the  elastic  restoring  force  is, 

FR  "  "  Chh 


C190) 
(191) 

(192) 

(193) 


the  sum  of  the  moments  of  the  inertia  forces  about  the 
elastic  axis  is, 

Mj.  =  -  (SQh  +  IQ0)  (194) 

and  the  elastic  restoring  moment  is, 

MR  =  -  C99  •  (195) 

L  and  Mfl  are  the  aerodynamic  lift  and  pitching  moment  about 

the  elastic  axis.   With  the  assumption  of  simple  harmonic 

oscillation,  the  blade  motion  can  be  written, 

h  =  h  e 
o 

.       z    ,       io)t 
h  =  -  0)   h  e 


6  =  0   e 
0 


o 

icot 


a  2   q     lWt 

e  =  -w    e_  e 


The  equations  of  motion  thus  become, 

L  =  eiwt[-  w2Mh   -  w2SQe   +  C,  h  J 

o  0o  ho 

M   =    eXWt[-    o)2S0h    -    a)2Io0      +    Cft6    ] 

o  *  00  0O  0     0 


(196a) 
(196b) 
(197a) 
(197b) 

(198) 
C199) 
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Equating  these  to  Garrick-Rubinow f s  standard  lift  and  moment 

expressions  and  following  Garrick-Rubinow '  s  (.1946)  method, 

the  equations  can  be  written  as, 

2h 

(fl.  X-y+  Ln+  i  L0) — 2.  +  (-yX 0+  LQ+  i  L,)0   =0  (200a) 

n       jl     z.      a.  yj     ho 

c 

2h 

(-yxQ  +  Mx  +  i  M2)— -  +(fiQX  -  yr2  +  M3  +  i  M4)0q=  0    (200b) 

c 

The  two  simultaneous  flutter  equations  form  an  eigenvalue 
2h 


problem  in 


and  0  ,  thus,  the  determinant  of  coefficients 
o 


(the  characteristic  determinant)  must  be  zero: 

-  yxQ  +  L-  +  i  L, 


fl.X  -  y  +  L,  +  i  L„ 
n  x      '  2. 


-yx  +  M   +  i  M2 


Sx"  yre  +  M3  +  i  M4 


=  0     (201) 


This  is  the  flutter  determinant  with  unknown  X.   If  the 
determinant  is  expanded,  and  the  real  and  imaginary  parts 
equated  to  zero  the  result  is  two  equations  that  given  a 
value  of  k  and  M  will  yield  values  of  X.   For  a  specific 
Mach  number,  the  flutter  point  occurs  at  the  value  of  k  that 
produces  identical  values  of  X  in  the  two  equations.   From 
this  value  of  k  and  X  the  flutter  speed  and  frequency  can  be 
calculated . 

The  two  flutter  equations  are: 
(Real  part) 

nhfiex2+[n0a1-u)  +  nh(M3-ur2)jx  t  cR  =  o  (202) 

(Imaginary  part) 

(neL2  +  \M4)X  +  Ci  =  °  t203) 

where, 
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CR   =    pIxeCWx   +   L3)-CM3-    yr*)~    Ljxl    -    VxzQ]f   DR  C204a) 


CI    =    Plx9CM2    +   L^   M4    -    L^JjfDj 


C204b) 


and   where, 


DR   «    LXM3    -    L3M1    -    L2M4    +    J.  fa  C205a) 

Dj    =    LXM4    -    hfa   +    L2M3    -    L^  (205b) 

Solving    for    X, 
(Real    Part) 

-[^e(L1-y)+^h(M3-yr^)]+/[^e(L1-y)+fih(M3-yr^]"-4CR^h^e 

2fih^6  (206) 


X   = 


(Imaginary  Part) 

X  =   "a  L?  I    57  M, (207) 

0   2     h   4 

The  solution  to  the  flutter  problem  is  then  to  obtain 
values  of  x  as  a  function  of  k  until  the  real  and  imaginary 
solutions  converge.   At  this  "flutter  point'  the  non-dimen- 
sional flutter  frequency  is, 

'Hi    jl. 

»6  "  /J"  (208) 

and  the  non-dimensional  freestream  flutter  speed  is 

(209) 


ctoQ   k/x 

It  should  be  noted  that  the  interblade  phase  angle  is  an 
open  parameter  (.as  well  as  k)  in  calculating  lift  and  moment 
on  the  cascade  blade.   Lane  (1956)  showed  that  the  minimum 
flutter  speed  of  a  cascade  occurs  at  a  non-zero  interblade 
phase  angle  and  that  this  is  the  critical  phase  angle  at 
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which  the  cascade  will  oscillate.   For  an  unstaggered  cascade 
this  critical  phase  angle  was  shown  to  be  180Q,  but  for  a 
staggered  cascade  it  can  only  be  found  by  repeated  flutter 
calculations  at  various  phase  angles.   The  problem  is  solved 
when  the  minimum  flutter  speed  is  determined. 

If  only  a  single  degree  of  freedom  of  motion  (pitch  os- 
cillation) exists,  the  flutter  problem  is  greatly  simplified. 
The  equation  of  motion  of  the  blade  can  then  be  written, 

iee  +  (i+ig)c0e  =  m0  (210) 

Here  the  sum  of  the  moments  of  the  inertia  forces  about  the 
elastic  axis  is, 

Mj  -  -  iQe  (2ii) 

the  elastic  restoring  moment  is, 

MR  =  -  (l+ig)CQ9  (212) 

and  the  aerodynamic  pitching  moment  about  the  elastic  axis 
is, 

MQ  =  -  7  p   c2u2k26  eiU)t(M0+  i  M.)  (213) 

9      4   °°    *»    o       3      4 

With  the  assumption  of  simple  harmonic  oscillation  as  before, 

the  blade  motion  can  be  written  as, 

G  =  6   eia3t  (214) 

o 

•  •  • 

9  =  -co2  6  elWt  (215) 

o 

The  equation  of  motion  then  becomes, 

eiWVl0u2eo+(l+ig)Ceep>  -  £pTOe2u2k2epeiWtCM3+i  M4)  C216) 

Separating  the  equation  into  its  real  and  imaginary  parts 
gives , 
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"V2    +   C9   +  JP^"^^   =    ° 
g   Ce   t  K2  VM4    "    Q 
or    as    in    Garrick    and    Rubinow    (1946) 
fi0X      -    Vr2Q    +   M3    -    0 


M4  +  g  ftQX  =  0 


(217a) 
(217b) 

(218a) 
(218b) 


Thus,  the  solution  to  the  single-degree-of -freedom  prob- 
lem is  to  obtain  values  of  M.  as  a  function  of  k  until  the 

4 

second  equation  is  satisfied.  The  value  of  M«  for  that  k  is 
then  used  in  the  first  equation  to  determine  if  X  is  reason- 
able : 


X  =  1  - 


M. 


>  0 


(219) 


a 


The  non-dimensional  flutter  frequency  is  as  before, 


0), 


we      /x 

and  the  non-dimensional  freestream  flutter  speed  is, 

U 


(220) 


_F 

CO) 


e 


k/x 


(221) 


Again  it  must  be  noted  that  in  the  case  of  staggered 
cascades  the  flutter  computations  must  be  repeated  with  vary- 
ing interblade  phase  angle  until  the  minimum  flutter  speed 
is  found.   This  will  give,  the  critical  flutter  condition. 

As  in  the  analysis  of  the  single  airfoil,  in  the  absence 
of  structural  damping,  aerodynamic . ins tability  (flutter) 
occurs  when  M,  ^  0.   When  M,  is  written  as  the  pitch  damping 
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derivative,  Eq.  (171b) ,  it  should  be  noted  that  aerodynamic 
undamping  occurs  when  the  pitch  damping  derivative  is  posi- 
tive. 
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VI.   COMPUTATIONAL  PROCEDURE 

The  computer  programs  used  to  calculate  the  cascade  flow 
field  are  based  on  one  used  by  Platzer  and  Pierce  (1970)  to 
calculate  pressure  distributions  on  an  airfoil  oscillating 
with  wind  tunnel  interference.   The  programs  here,  however, 
in  addition  to  calculating  the  pressure  distribution  on  the 
blade  of  a  cascade  integrate  the  pressure  distribution  over 
the  blade  surface  and  calculate  the  flutter  frequency  and 
freestream  flutter  velocity.   One  program  (Program  A)  is  for 
two-degree-of-f reedom  (pitch  and  plunge)  flutter  calculations 
while  the  second  (Program  B) ,  merely  a  simplification  of  the 
first,  is  for  single-degree-of -freedom  (pitch)  flutter  cal- 
culations . 

Both  programs  are  valid  only  for  cascades  with  supersonic 
leading-edge  locus  vibrating  with  small  amplitude  oscilla- 
tions.  The  programs  calculate  the  flow  field  between  two 
adjacent  blades  using  the  method  of  characteristics  finite 
difference  equations  developed  previously  and  integrate  the 
pressure  distributions  on  the  blade  surfaces  using  the  trape- 
zoidal rule. 

The  computational  molecules  shown  in  Figures  2,  3,  and  A 
are  used  to  calculate  the  flow  field  quantities  at  P„~  in 
the  general  flow  field,  on  the  lower  blade,  and  on  the  upper 
blade,  respectively.   The  distance  Ax  shown  in  the  figures  is 
determined  by, 
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Ax  ^  ~  '  (222) 

where  v  is  the  grid  fineness  ratio.   This  parameter  is  an 
arbitrary  input  variable  into  the  programs  and  is  equal  to 
one  less  than  the  total  number  of  grid  points  on  any  Mach 
line  that  runs  from  one  blade  to  the  other. 

Input  parameters  are  entered  into  the  programs  by  means 
of  the  NAMELIST  option  of  FORTRAN.   In  using  this  option,  in- 
put parameters  are  combined  in  a  titled  list  and  then  referred 
to  in  the  programs  by  this  title.   In  both  Program  A  and  B 
the  NAMELIST  title  is  the  same  (NAM1) ;  however,  the  input 
variables  in  each  are  slightly  different.   The  real  advantage 
of  using  NAMELIST  is  that  each  data  card  need  only  contain 
the  value  of  one  input  variable  while  the  format  of  the  data 
card  is  such  that  the  variable's  name  is  punched  on  the  card 
as  well  as  its  value.   In  this  way,  when  a  different  value  for 
an  input  variable  is  desired,  rather  than  changing  the  entire 
data  deck,  only  one  easily  identifiable  data  card  need  be 
changed . 

The  format  for  a  NAMELIST  data  deck  is  as  follows:   The 
first  column  in  each  card  of  the  deck  is  left  blank.   The 
first  card  starts  in  the  second  column  with  the  symbol  & 
(ampersand)  followed  immediately  by  the  title  of  the  NAMELIST 
(here,  NAM1) .   The  next  cards  list  the  variables  and  their 
input  values,  one  per  card,  in  any  order.   The  format  here, 
starting  in  the  second  column  and  followed  by  a  comma,  is: 
Variable  Name  =  Value.   The  last  card  in  the  data  deck  and 
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the  signal  to  the  computer  that  the  NAMELIST  has  ended  again 
starts  in  the  second  column  with  the  symbol  &  Campersand) 
followed  immediately  by  END. 

The  value  of  an  input  variable  may  be  written  in  any  for- 
mat as  long  as  integer  variables  are  written  without  decimal 
points  and  real  variables  are  written  with  a  decimal  point. 
Blanks  are  taken  as  zeroes,  but  a  comma  must  appear  somewhere 
between  a  desired  input  value  and  the  next  data  card. 

In  both  programs,  prior  to  the  NAMELIST  dataset,  the  date 
is  entered  on  a  data  card  in  the  first  twelve  spaces. 

To  illustrate  the  input  procedures,  the  following  two 
examples  are  proposed: 

Two  Degree  of  Freedom  Flutter   One  Degree  of  Freedom  Flutter 
v  =  100  v  =  100 

M  =  /I  M  =  /I 

Y  =  1.4  y    =  1.4 

k  =  2  k  =  2 

xQ  =  0.5  xQ  -  0.5 

d  =  0.6  d  =  0.6 

3  =  26.565°  3  =  26.565° 

6  =  180°  6  =  180° 

y  =  500  y  =  500 

r0  =  0.5  *0  "  °*5 

xA  =0.1  g  =  2% 


~T   "'1.5  A2/k  m    1 


0 

A2/k  =1  8=1.2 

c  =  1.2 
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The  data  cards 
Card   Program 

1  A  &  B 

2  A  &  B 

3  A  &  B 


needed  are  as  follows : 


4 

A  & 

B 

5 

A  & 

B 

6 

A  & 

B 

7 

A  & 

B 

8 

A  & 

B 

9 

A  & 

B 

10 

A  & 

B 

11 

A  & 

B 

12 

A  & 

B 

13 

A 

13 


B 


Format 
16  JUNE  1972 
&NAM1 
NGRDFN  =  100, 


Description 

Date 

NAMELIST  title 

Grid  fineness  ratio,  must 
be  an  even  integer,  less 
than  400. 


FSTRMN  =  1.414214,  Freestream  Mach  number 
RTOSPH  =  1.4,       Ratio  of  specific  heats 

Reduced  frequency 
Elastic  axis  position 
Blade  distance 


REDFRQ  =  2.0, 
XSUB0  =  0.5, 
TNWDST  =  0.6, 
STGANG  =  26.565, 


FAZE  =  180.0, 
MUU  =  500.0, 
RSUBA  =  0.5, 
XSUBA  =  0.1, 

GSUBA  =  0.02, 


14 

A  & 

B 

INCRE  =  1.0, 

15 

A  & 

B 

FIN  =  1.2, 

16 

A 

HAFREQ  =  1.5, 

16 

B 

&END 

17 

A 

&END 

Stagger  angle,  tan  3^ 
cot  a. 

Interblade  phase  angle 

Blade  density  parameter 

Radius  of  gyration 

Distance  from  the  blade 
elastic  axis  to  its  center 
of  gravity. 

Structural  damping  coeffi- 
cient 

2/k  increment 

Chord  length,  c  ^  2A. 

Ratio  of  natural  frequencies, 
bending  to  torsion. 

End  of  NAMELIST 

End  pf  NAMELIST 
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The  definitions  of  these  parameters  are  given  again  in  the 

alphabetical  listing  of  all  program  variables  in  Appendix  A. 

Program  A  has  eight  subroutines.   They  are  as  follows: 

INPUT     Reads  in  all  input  data  and  sets  up  the 
finite  difference  grid. 

INTIAL    Initializes  all  the  flow  quantities  at  the 
leading-edge,  starts  the  integration,  and 
initializes  most  of  the  logic  variables. 

MACHLN    Computes  the  flow  quantities  along  the 

initial  left  and  right-running  Mach  lines 
at  the  given  grid  point. 

HIFOIL  Computes  the  flow  quantities  on  the  lower 
surface  of  the  upper  airfoil  at  the  given 
grid  point,  and  continues  the  integration 
of  the  pressure  there. 

GENFPT  Computes  the  flow  quantities  at  a  general 
flow  field  grid  point. 

LOFOIL  Computes  the  flow  quantities  on  the  upper 
surface  of  the  lower  airfoil  at  the  given 
grid  point,  and  continues  the  integration 
of  the  pressure  there. 

COMPXY    Computes  the  x  and  y  position  of  the  given 
grid  point . 

FLUTER    Computes  the  oscillation  frequency  ratio, 
X,  for  both  real  and  imaginary  flutter 
equations . 

Program  B  contains  all  of  these  subroutines  except  FLUTER. 
Since  the  flutter  problem  is  greatly  simplified  in  the  single 
degree  of  freedom  case,  the  flutter  calculations  are  made  in 
the  main  program. 

Program  A  contains  at  least  two  "DO"  loops.   The  inner 
one  is  for  calculating  the  flow  field  with  first  pitch  oscil- 
lation boundary  conditions  and  then  plunge  oscillation 
boundary  conditions.   This  loop  is  not  present  in  Program  B 
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since  only  the  flow  field  due  to  pitch  oscillation  is  desired. 
The  next  outer  loop  is  used  to  increment  2/k  (1/k  in  Garrick- 
Rubinow's  (1946)  notation)  starting  at  that  reduced  frequency 
specified  in  INPUT.   This  loop  is  continued  until  hopefully 
the  flutter  point  is  found.   Any  loops  outside  of  these  are 
used  to  make  multiple  flutter  calculations  in  one  run  with 
varying  input  parameters. 

In  INPUT  the  grid  is  determined  using  the  grid  fineness 
ratio,  blade  spacing,  and  Mach  number.   The  value  of  Ax  as 
shown  in  Eqs.  (96)  is  obtained.   The  stagger    angle  is  then 
adjusted  such  that  the  distance  the  upper  blade  is  staggered 
back  is  the  nearest  multiple  number  of  Ax  increments.   The 
"Compatible  Stagger  Angle"  is  then  printed  out.   The  upper 
blade  may  be  staggered  back  such  that  the  distance  between 
its  leading-edge  and  its  intersection  with  the  initial  left- 
running  Mach  line  from  the  lower  blade  is  zero,  but  no  further 
With  less  stagger  this  distance  must  be  greater  than  or  equal 
to  2Ax,  or, 

A  -  B  £  2Ax  (223) 

This  allows  a  minimum  of  three  grid  points  in  the  primed  flow 
field  zones  as  shown  in  Figure  5.   This  is  the  minimum  number 
of  points  required  when  the  value  of  the  last  point  must  be 
linearly  extrapolated  from  the  other  two,  as  is  the  procedure 
in  the  programs. 

Flow  quantities  along  the  initial  left  and  right-running 
Mach  lines  are  computed  using  the  equations  derived  in  the 
Method  of  Characteristics  for  initial  conditions.   Those 
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quantities  along  the  blade  surfaces  are  computed  using  the 
equations  derived  for  the  appropriate  boundary  condition. 
The  values  of  the  normal  velocities  are  derived  as  a  result 
of  the  analysis  of  the  flow  tangency  condition  and  have  been 
given  previously. 

Because  the  programs  compute  the  flow  field  in  rows  of 
right-running  Mach  lines  starting  with  points  on  the  initial 
left-running  Mach  line,  provision  must  be  made  for  the  com- 
putation of  points  in  the  area  referred  to  as  Zone  I'.   The 
reason  is  that  once  the  computation  reaches  the  intersection 
of  the  initial  left-running  Mach  line  from  the  lower  blade 
and  the  initial  right-running  Mach  line  from  the  upper  blade, 
the  points  on  the  initial  left-running  Mach  line  become 
general  flow  field  points  requiring  those  equations  for  the 
flow  quantity  computation.   To  use  those  equations  the  flow 
quantities  along  the  last  left-running  Mach  line  in  Zone  I' 
must  be  known.   They  are  obtained  in  two  ways  depending  on 
the  value  of  the  interblade  phase  angle:   If  the  phase  angle 
is  non-zero,  the  entire  flow  field  in  Zone  I'  (along  with  the 
integration  over  the  bottom  surface  of  the  top  blade)  is  com- 
puted concurrently  with  the  computation  of  the  flow  field  in 
Zone  I.   If  the  phase  angle  is  zero,  only  that  row  in  Zone  I 
that  is  the  image  of  the  last  left-running  Mach  line  in  Zone 
I'  must  be  determined  since  the  flow  quantities  of  these  two 
rows  are  the  same  (.except  for  the  signs  of  U  and  C)  . 

The  pressure  distributions  on  the.  surfaces  of  the  blades 
are  integrated  point  by  point  as  they  are  being  calculate 
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Since  the  pressure  coefficient  at  the  last  point  in  each  zone 
of  influence  is  unknown,  its  yalue  must  be  linearly  extrapo- 
lated from  the  preceding  points.   The  pressure  distribution 
of  the  upper  blade  is  adjusted  so  as  to  represent  the  pressure 
distribution  on  the  bottom  surface  of  the  lower  blade: 

C(x,0~)  =  e"i<5  C(x+B,d)  (224) 

In  this  way  the  result  of  the  integration  is  the  lift  and 
pitching  moment  on  a  representative  blade  of  the  cascade. 

Termination  of  the  flow  field  calculation  starts  when  the 
distance  from  the  leading-edge  of  the  lower  blade  to  the 
point  of  calculation  on  the  lower  blade  reaches  the  prescribed 
chord  length.   Because  HIFOIL  must  be  called  at  least  once  for 
proper  program  termination,  this  chord  length  has  a  lower 
bound : 

c  ^  2A  (225) 

Since  the  flow  field  is  being  calculated  along  rows  of  right- 
running  Mach  lines,  the  last  zone  along  the  upper  blade  is 
then  computed  to  complete  the  calculation.   Although  the  flow 
field  is  terminated  when  the  calculation  reaches  the  prescribed 
chord  length,  the  integration  of  the  pressure  distributions 
is  terminated  when  the  chord  length  reaches  unity. 

When  the  values  of  total  lift  and  pitching  moment  have 
been  determined,  Program  A  enters  them  in  the  flutter  equa- 
tions, Eq.  (206)  and  Eq.  C2Q7)  via  subroutine  FLUTER.  Program 
B,  computing  only  the  pitching  moment,  enters  these  values  in 
Eqs.  (218).   Convergence  to  the  flutter  point  is  then  deter- 
mined . 
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In  Program  A,  the  difference  between  solutions  to  the 
real  flutter  equation  and,  that  to  the  imaginary  flutter  equa- 
tion is  noted  after  each  computation.   If  the  absolute  value 
of  the  difference  is  less  than  10   ,  the  flutter  point  is 
considered  to  be  determined.   When  the  difference  changes  sign, 
indicating  convergence  somewhere  in  the  frequency  interval, 
the  2/k  increment  is  divided  by  five  and  the  computation  is 
started  again  at  a  new  value  of  2/k  determined  by  adding  the 
new  increment  to  the  previous  lower  value.   Each  time  the 
flutter  point  is  passed  the  increment  is  made  smaller.   The 
flutter  point  is  considered  to  be  ascertained  when  the  incre- 
ment is  less  than  or  equal  to  10 

In  Program  B,  the  value  of  X  is  determined  from  Eq .  (218a). 
This  value  is  then  used  along  with  M,  to  determine  the  value 
of  Eq.  (218b).   If  the  absolute  value  of  Eq .  (218b)  is  less 
then  10   ,  the  flutter  point  is  considered  to  be  determined. 
If  Eq.  (218b)  is  negative,  indicating  convergence  somewhere 
in  the  frequency  interval,  the  same  procedure  is  followed  as 
in  Program  A.   When  the  flutter  point  is  determined,  the 
value  of  X  is  tested  to  insure  a  realistic  flutter  frequency 
as  shown  in  Eq.  (219). 

The  output  from  Program  A  consists  of  the  input  parameters 
(including  such  items  as  x  and  the  compatible  stagger  angle); 
the  values  for  X   and  X   as  defined  in  Eq.  C206)  and  Eq .  (207), 
respectively,  for  each  reduced  frequency  calculated;  and  the 
non-dimensional  flutter  frequency  and  flutter  speed  as  defined 
in  Eq.  (208)  and  Eq .  (209),  respectively.   The  value  of 
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"1/k"  given  in  the  output  corresponds  to  the  reciprocal  of 
the  reduced  frequency  as  defined  by  Garrick  and  Rubinow 
(1946).   This  is  equivalent  to  our  2/k. 

The  output  from  Program  B  is  identical  to  that  of  Pro- 
gram A  except  that  Program  B  gives  the  values  of  M , ,  M  ,  X, 
and  the  system  damping  (M,  -  gfi„X)  as  shown  in  Eq .  (218b) 
for  each  reduced  frequency  calculated. 

A  flow  diagram  of  Program  A  is  given  in  Appendix  B. 
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yil.   RESULTS 

A.   COMPARISON  OF  ELEMENTARY  THEORY  WITH  MILES 

Miles  (1956)  presented  an  analysis  of  supersonic  flow 
past  an  oscillating  airfoil  subjected  to  wall  interference 
in  a  solid-wall  wind  tunnel.   Using  Laplace  transform  tech- 
niques he  derived  the  following  results  for  the  moment  co- 
efficient (his  Eq.  (7.18)) 

Cmct=  4  tan  a{(2N+l)[x  -  ^]-N(N+l)Ax  +  ^N  (N+l)  (2N+1)  A2 }  (226a) 

O    Z  O    o 

c   =  -  4  tan  a{(2N+l)[x2-  x  +  \\    -    N(N+1)A 
mq  o    o   J 

[x2-  x  +  h    +   ^-N2(N+1)2A3}  (226b) 

o    o   z     lz 

c  .=  4  tan  a{-(2N+l)  tan2a[^  x   -  ^-]-N(N+l)A 
ma  2   o    3 

[x  -  ~]+  \   N(N+1) (2N+1)  (tan2a  +2) 

O    Z      D 

•  A2xq-  ^yN2(N+'l)2(2  tan2a  +  3)A3}       (226c) 

where  N  is  the  largest  integer  smaller  than  A    and  a,  d 

subscripts  represent  angle  of  attack.   Here,  for  a  slowly 

oscillating  airfoil, 

C   =  a  [c   +  ik(c    +  c  .)]e±kt  (227) 

m     o   ma       mq     ma 

unlike  Eq.  (170)  . 

As  stated  by  Miles,  this  case  is  identical  to  that  of 

the  unstaggered  cascade  oscillating  with  an  interblade  phase 

angle  of  180°,  and  the  yalues  of  pitching  moment  as  defined 

in  Eq.  (.170)  and  Eq .  (.171)  based  on  our  elementary  theory 

for  B  -  0,  5  -  180,  and  0.5  <  A  ^  1.0-,  prove  to  be  the  same 

as  those  shown  in  Eqs .  (226)  with  N  =  1.   The  torsional 
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stability  boundaries  c  *  m    0  for  an  unstaggered  cascade  os- 

*  Bit? 

dilating  with  an  interblade  phase  angle  of  180°  as  computed 
from  Eq.  (171b)  are  shown  in  Figure  6  for  various  values  of 
solidity. 


A  =  0.7 


A  =  0.9 


Unstable  (A  <  l) 


A  =  1.0 


Unstable 

I 1 I L 


.1    .2    .3    .k         .5    .6    .7  Xq 
Figure  6 

The  large  destabilizing  influence  of  interference  causing 
instability  over  the  whole  supersonic  Mach  number  range  for 
a  wide  range  of  elastic  axis  positions  can  be  seen.   Also 
shown  is  the  interference  -  free  boundary  (A  =  1)  for  the 
single  airfoil  in  an  unbounded  supersonic  flow  (Garrick  and 
Rubinow  1946)  exhibiting  instability  as  is  well  known  only 
for  low  supersonic  Mach  numbers.   As  demonstrated  above, 
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these  stability  boundaries  also  apply  to  the  slowly  oscillat- 
ing airfoil  subjected  to  solid-wall  wind  tunnel  interference. 

B.  COMPARISON  OF  ELEMENTARY  THEORY  WITH  LANE 

Lane  (1957)  extended  Miles'  Laplace  transform  solution  to 
the  analysis  of  oscillating  cascades  with  supersonic  leading- 
edge  locus.   By  expanding  the  blade  pressure  distributions 
(his  Eq .  (13))  with  respect  to  frequency  and  retaining  only 
terms  up  to  the  first  power  in  frequency,  the  pressure  jump 
across  a  slowly  oscillating  blade  can  be  obtained.   This 
pressure  jump  is  identical  to  that  being  integrated  in  Eq . 
(168)  . 

C.  COMPARISON  OF  ELEMENTARY  THEORY  WITH  DRAKE 

Drake  (1956)  using  Laplace  transform  techniques  presented 

an  analysis  of  supersonic  flow  past  an  oscillating  airfoil  in 

a  wind  tunnel  with  free  jet  boundaries.   His  pitching  moment 

coefficient  is  defined  in  the  English  system,  equivalent  to 

Eq.  (170),  is, 

C   =  a  [m   +  ik  m-]  elkt  (228) 

m     o   a       a  v    ' 

for  a  slowly  oscillating  airfoil  where  a  again  refers  to 

angle  of  attack.   Here, 

ma    =    4    tan    aUxQ-l)h1CA)+    h2(A)J  C229a) 


m.    «    4    tan   a    {h~(A)-   h0(A)-x    (x   -l)h, (A) 

u  O  Z  Q         0  JL 


f    tanza    (M2A  f^   -    1) [ {xQ-l}h2 (A)t   h3(A)J}  C229b) 

where    for    0.5   ^   A   ^  1.0, 
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o 
mome 


h^U)  *  2A  -  1  C230a) 

K2CA)  m    -A2  +  2A  -  \  C230b) 

h3CA)  =  jA3  -  A2  +  A  -  |  (230c) 

This  case  is  identical  to  that  of  the  unstaggered  cascade 
scillating  in  phase  (6=  0) ,  and  the  values  of  pitching 

nt  as  defined  in  Eq.  (170)  and  Eq.  (171)  based  on  our 
elementary  theory  for  6=0,  3=0,  and  0.5  <  A  ^  1.0  prove 
to  be  the  same  as  those  shown  in  Eq.  (229). 

Drake  (1957)  extended  his  work  to  the  case  of  the  airfoil 

oscillating  in  a  wind  tunnel  with  porous  walls.   Platzer 

(1971)  will  publish  an  analysis  of  this  problem  using  our 

elementary  theory.   His  pitch  damping  coefficient  as  defined 

in  Eq.  (170),  for  a  slowly  oscillating  airfoil  is, 

c  *  =  -  4  tan  a{x2  +  x  (■=■  tan2a-l)+  t(1-  tan2a) 
mo  o     o  L  5 

2(1-A)   j-  r ..  „       Wn     2  n.  2a  tan3a, 

+  tt — i — —  I  [1-0  tan  a   1-tan  a  +  r— J 

1+atana   l  1+Ctan  a 

[X]  +  [(1-0    tan  a) (xq+  A(l+tan2a) ) ] [ Y ] } }       (231) 

where 

X  =  (1-A) [-  jxq+   |(A+2] 

Y  =  xq  -  |  (A+l) 

and  0.5  <  A  ^  l.Q.   Here  Q  <  a  <  "°  where  a  =  0  represents  the 
solid-wall  wind  tunnel  case  and  a  =  °°  represents  a  free-jet 
wind  tunnel  boundary.   With  these  limiting  values  of  porosity 
substituted  in  Eq .  (231)  the  comparison  with  Eq.  (171b)  is 
exact  for  the  case  of  6  -  180  and  (5  =  0,  respectively. 
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Further  the  results  check  with  Drake  Q-957).   Stability  boun- 
daries for  an  airfoil  pscillating  slowly  in  a  porous-wall 
wind  tunnel  are  given  in  Figure  7  for  a  tunnel  aspect  ratio, 
A  =  0.5  and  for  various  values  of  tunnel  porosity.   An  in- 
crease in  tunnel  porosity  is  seen  to  have  a  large  stabilizing 
influence . 


D.   AIRFOIL  -  AIRFOIL  INTERFERENCE 

Our  elementary  solution  pf  the  slowly  oscillating  cascade 
problem  can  be  used  to  study  pther  interference  problems  as 
well.   The  interference  between  two  airfoils  may  be  of 
interest  for  certain  proposed  space  shuttle  configurations, 
although  two-dimensional  flow  can  only  give  an  indication  of 
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trends.   Our  elementary  solution  applied  to  the  unstaggered 
cascade  oscillating  with  an  interblade  phase  angle  of  180°, 
as  stated  before,  is  equivalent  to  an  airfoil  subjected  to 
solid-wall  wind  tunnel  interference.   If  the  upper  wind  tun- 
nel wall  is  considered  a  large  stationary  airfoil  (idealized 
as  a  flat  plate  at  zero  angle  of  attack),  the  theory  can  be 
used  to  determine  the  non-dimensional  pressure  distribution 
along  the  upper  surface  of  a  smaller  airfoil  mounted  closely 
below  it.   Using  the  non-interference  pressure  distribution 
over  the  lower  surface  of  the  airfoil,  the  net  pressure  dis- 
tribution over  the  airfoil  can  be  integrated  over  the  chord 
length  and  the  pitch  damping  coefficient  as  defined  in  Eq . 
(171b)  can  be  determined  as, 

c  •  =  -2  tan  a  {x2(4-2A)-  x  [3A2-  4A+  2 
m6  o  o 

+  (A2-  2)(tan2a-l)]  +  2A3-2A  -  |(2-A3) 

(tan2  a  -1)}  (232) 

where 

A  =  2h  cot  a 
h  =  airfoil  distance 
and  0.5  ^  A  <  1.0. 

The  stability  boundaries  for  an  airfoil  undergoing  inter- 
ference of  this  type  are  shown  in  Figure  8.   The  effect  of 
interference  is  strongly  destabilizing.   This  case  further 
applies  to  an  airfoil  mounted  close  to  only  one  solid  wall 
in  a  wind  tunnel. 
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E.   COMPARISON  OF  THE  METHOD  OF  CHARACTERISTICS  WITH 
ELEMENTARY  THEORY 

Values  of  flow  field  quantities,  U,  V,  and  C,  both  real 
and  imaginary,  based  on  computations  using  Program  B  were 
compared  with  those  quantities  using  the  elementary  solution 
as  applied  to  cascades.   Comparisons  were  made  for  various 
values  of  A,  stagger,  and  interblade  phase  angle.   Points 
along  both  surfaces  of  the  blades  as  well  as  in  the  flow 
field  were  compared.   In  the  case  of  zero  interblade  phase 
angle,  values  of  A  as  low  as  0.25  were  compared.   The  restric- 
tion of  supersonic  leading-edge  locus  was  maintained.   All 


75 


values  compared  were  wi.th.in  2%  up  to  reduced  frequencies  of 
0.1.   Differences  between  1%  and  2%  were  obtained  only  at 


Mach  numbers  near  one  CM  w  /l . 2 5 )  at  a  reduced  frequency  of 
0.1.   A  sample  of  the  flow  quantities  compared  is  given  in 
Table  1. 


C22R 

C22R 

C22I 

C22I 

// 

X 

y 

(M.O.C.) 

(S  .S.C.) 

(M.O.C.) 

(S.S.C.) 

0.10 

0.00 

-1.000 

-1.000 

0.000 

0.000 

0.50 

0.00 

0.659 

0.660 

1.120 

1.120 

1 

0.60 

0.00 

0.659 

0.660 

1.120 

1.120 

0.10 

0.50 

0.800 

0.800 

0.600 

0.600 

0.50 

0.50 

-1.197 

-1.200 

0.699 

0.700 

0.60 

0.50 

-1.197 

-1.200 

0.699 

0.700 

0.20 

0.00 

-1.000 

-1.000 

0.000 

0.000 

0.50 

0.00 

0.659 

0.660 

1.119 

1.120 

2 

0.60 

0.00 

0.659 

0.660 

1.119 

1.120 

0.44 

0.50 

0.800 

0.800 

0.600 

0.600 

0.74 

0.50 

-1.197 

-1.200 

0.699 

0.700 

0.84 

0.50 

-1.197 

-1.200 

0.699 

0.700 

0.20 

0.00 

-1.997 

-2.000 

0.060 

0.060 

0.50 

0.00 

-1.984 

-2.000 

0.149 

0.150 

3 

0.74 

0.00 

1.510 

1.500 

2  .203 

2.222 

0.80 

0.00 

1.530 

1.522 

2.187 

2.211' 

0.40 

1.00 

1.627 

1.629 

1.161 

1.162 

0.50 

1.00 

-2.322 

-2.353 

1.633 

1.650 

0.74 

1.00 

-2.252 

-2.264 

1.708 

1.738 

in 

M    = 

/2 

n 

M    =    /T 

#3 

M   =    /l. 25 

A    = 

0.5 

A    =    0.5 

A   =    0.5 

3   = 

0.0 

3    =    25.6 

4 

3    =    13.50 

<5   - 

36.87 

<$    n    36.8 

7 

6    -    36.87 

x    - 

0 

0.0 

X    ^    0.0 

0 

x    *    0.0 

o 

Table  1 
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Further,  comparisons  of  the  out  of  phase  pressure  distribu- 
tion on  the  upper  surface  of  a  staggered  cascade  blade  in 
supersonic  flow  of  M  -  /l . 25  and  M  f  v2    as  predicted  by  the 
method  of  characteristics  and  the  elementary  theory  are 
shown  in  Figure  9  (M  -  /2,  B  =  .12,  xQ  -  0,  6  =  0)  and 


Figure  10  (M  =  /1.25,  B  =  .12,  x   =  0,  5  =  0).   Here,  the 
solid  lines  indicate  the  elementary  theory,  and  the  dashed 
lines  indicate  the  method  of  characteristics. 
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Figure  9 

Agreement  between  the  two  methods  is  quite  good  for  re- 
duced frequencies  up  to  0.1  at  H  ■  /2~,  but  this  close  agree- 
ment is  maintained  at  M  =  /l . 25  only  up  to  a  reduced  frequency 
of  0.05.   Since  both  methods  are  based  on  the  same  linearized 
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equations  of  motion,  this  shows  that  the  influence  of  fre- 
quency is  greater  at  the  iower  supersonic  Mach  numbers. 

Integrated  values  of  lift  and  moment  (Xo  and  M.)  were 
also  compared.   Here  differences  as  high  as  3%  were  found  at 
Mach  numbers  approaching  one;  however,  for  M  =  /l ,  the 
values  were  within  2%. 


i  [c  /e  ] 

m  p  o J 


Figure  10 

These  results  were  produced  using  a  grid  fineness  ratio 

of  100  or  one  quarter  of  the  capacity  of  the  program  in  order 

to  reduce  the  amount  of  computer  time  required  to  make  the 
computations . 
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F.  COMPARISON  OF  METHOD  OF  CHARACTERISTICS  WITH 
GARRICK-RUBINQW 

Lift  and  moment  calculations  obtained  from  Programs  A 
and  B  at  A  =^  1  were  compared  with  those  of  Garrick  and  Rubi- 
now  (1946)  for  both  pitch  and  plunge  oscillation.   Differ- 
ences were  less  than  2%.   Program  B  was  used  to  duplicate 
several  different  single-degree-of -freedom  flutter  points 
shown  in  Garrick  and  Rubinow's  (1946)  Figure  22.   Again,  the 
grid  fineness  ratio  used  was  100. 

G.  METHOD  OF  CHARACTERISTICS  RESULTS 

The  following  results  were  obtained  using  Program  B  to 
compute  single-degree-of -freedom  flutter  speeds  and  frequen- 
cies.  Approximately  four  minutes  of  CPU  time  on  the  IBM-360 
were  required  to  generate  one  flutter  point.   A  grid  fineness 
ratio  of  100  was  used  throughout. 
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Figure    11 
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In  Figure  11  CM  =  1.11,  A  =  Q.95,  <5  =  0,  xQ  =  .5,  u  =  15.708, 
Xq    =  .5,  g  ■  0)  a  typical  variation  of  M,  with  1/k  is  shown 
for  an  unstaggered  cascade  oscillating  in  phase  with  only  a 
slight  amount  of  blade  interference.   Here,  the  flutter  point 
is  at  a  relatively  low  value  of  reduced  frequency  (k  =  0.02). 
That  this  is  not  the  critical  flutter  condition  is  shown  in 
Figure  12  where  the  non-dimensional  flutter  speed  is  plotted 
versus  interblade  phase  angle  for  the  same  parameters  as  in 
Figure  11.   In  complete  agreement  with  Lane  (1956),  the 
critical  flutter  condition,  i.e.,  the  lowest  flutter  speed, 
for  this  unstaggered  cascade  is  at  an  interblade  phase  angle 
of  180°. 
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Figure  13  shows  the  effect  of  a  small  amount  of  interfer- 
ence on  the  flutter  speed  at  M  =  1.111  when  plotted  against 
elastic  axis  position. 
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Figure  13 

Again  the  blade  and  flow  parameters  are  those  as  shown  in 
Figure  11  and  Figure  12.   It  can  be  seen  that  a  relatively 
small  amount  of  interference  significantly  lowers  the  non- 
dimensional  flutter  speed  over  the  entire  range  of  elastic 
axis  positions  tested. 

In  Figure  14  (A  =  .95,  6  -  180,  B  -  0,  u  =  500,  g  =  0, 
rfl  =  .5)  the  effect  of  Mach  number  on  flutter  speed  is  shown 
for  various  elastic  axis  positions.   An  increase  in  Mach  num- 
ber to  M  =  /l    significantly  lowers  the  non-dimensional  flutter 
speed  for  a  relatively  small  amount  of  interference. 
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In  Figure  15  (6  =  180,  y  =  15.708,  rQ  =  .5,  g  =  0)  the 
effect  of  this  relatively  small  amount  of  interference  on 
the  flutter  speed  is  again  shown,  but  here  at  the  higher 
Mach  number  (M  =  v2 ) .   As  before,  the  non-dimensional  flut- 
ter speed  is  significantly  lowered  over  the  entire  elastic 
axis  range  by  a  small  amount  of  interference. 
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Figure  16  (M  =  /2,  A  =  .95,  B  *    0,  «5  =  180,  x   =  .5, 

o 

=  .5)  shows  the  effect  of  structural  damping  on  the  non- 


dimensional  flutter  speed  for  two  values  of  wing  density 
parameter . 
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Figure  16 

With  a  relatively  small  amount  of  interference,  the  flutter 
speed  is  greatly  increased  by  a  large  wing  density  parameter 
and  high  structural  damping.   The  high  value  of  wing  density 
parameter  is  more  applicable  to  the  cascade  case.   As  defined 
in  Garrick  and  Rubinow  (1946),  its  increase  may  be  interpreted 
as  an  increase  in  altitude  for  a  fixed  wing  density,  high 
values  indicating  supersonic  wings  at  high  altitude.   For 
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typical  supersonic  cascades,  the  wing  density  parameter  is 
of  the  order  of  700. 

The  remaining  results  give  variations  in  non-dimensional 
flutter  speed  of  the  staggered  cascade  ((3  -  26.565)  at  M=/2. 
Results  shown  are  based  on  a  cascade  with  an  interblade  phase 
angle  of  180°.   It  must  be  re-emphasized  that  this  interblade 
phase  angle  is  not  necessarily  the  one  that  gives  the  criti- 
cal flutter  condition;  however,  several  points  at  random  were 
chosen  and  the  critical  flutter  condition  found  for  each. 
All  were  within  1%  of  flutter  speed  at  5  =  180  and  within 
10°  of  6  =  180.   This  was  not  done  for  all  points  shown  in 
an  effort  to  reduce  the  total  amount  of  computer  time  used. 
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In  Figure  17  (y  =  500,  g  =  .02,  rQ  *  .5)  the  effect  of  an 
increasing  amount  of  interference  on  flutter  speed  is  shown 
for  various  elastic  axis  positions.   The  cascade  is  assumed 
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to  have  a  high  wing  density  parameter  and  2%  structural  damp 
ing  coefficient.   An  increase  in  the  amount  of  interference 
is  seen  to  lower  the  non-dimensional  flutter  speed. 
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Figure  18  (A  =  .6,  r„  =  .5)  shows  the  effect  on  the  non- 

u 

dimensional  flutter  speed  of  increasing  the  wing  density 
parameter  then  the  structural  damping  for  various  elastic 
axis  positions.   Without  structural  damping,  an  increase  in 
wing  density  parameter  had  little  effect  on  the  flutter 
speed.   With  2%  structural  damping,  however,  the  flutter 
speed  was  greatly  increased.   The  effect  of  structural  damp- 
ing is  shown  further  in  Figure  19  CA  -  .6,  x~  =  .5,  rD"  .5). 

u         v 

Here,  the  variation  of  flutter  speed  with  structural  damping 
is  shown  for  three  values  of  wing  density  parameter.   As 
can  be  seen,  increasing  the  structural  damping  increases  the 
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difference  in  flutter  speeds  for  the  three  values  of  wing 
density  parameter.   Thus,  the  accurate  determination  of  the 
wing  density  parameter  becomes  increasingly  more  important 
with  increased  structural  damping. 
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Figures  20  and  21  show  the  same  parameters  as  Figures  18 
and  19,  respectively,  but  here,  there  is  more  interference 
(A  =  .25).   The  results  are  relatively  the  same  as  those 
shown  in  Figures  18  and  19.   Flutter  speeds  overall  are  lower 
due  to  the  increase  in  the  amount  of  interference,  and  the 
effect  of  structural  damping  is  more  pronounced. 
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VIII.   CONCLUSIONS  AND  RECOMMENDATIONS 

The  elementary  theory  for  cascades  with  supersonic  lead- 
ing-edge locus  is  valid  for  low  frequencies  and  can  be  used 
to  give  rapid  calculations  of  lift  forces  and  pitching  mom- 
ents of  cascades  as  well  as  their  torsional  stability  boun- 
daries. These  results  can  immediately  be  applied  to  a  wide 
variety  of  wind  tunnel  and  airfoil  interference  problems. 

The  method  of  characteristics  can  be  used  in  conjunction 
with  the  high-speed  computer  to  give  rapid,  accurate  flutter 
calculations  of  cascades  with  supersonic  leading-edge  locus 
oscillating  at  arbitrary  frequency.   This  applies  to  torsion 
and  bending  flutter  calculations  as  well  as  to  torsion  flut- 
ter alone,  for  the  only  differences  in  the  calculations  are 
those  due  to  different  values  of  normal  velocity  at  the  blade 
surface.   In  each  case,  torsion  or  bending,  the  method  of 
calculating  the  flow  field  and  determining  the  lift  forces 
and  moments  is  identical. 

Results  of  single-degree-of -freedom  (torsional)  flutter 
calculations  show  that  interference  on  cascade  blades  due  to 
reflected  Mach  waves  from  adjacent  blades  lowers  flutter 
speeds  significantly  and  in  proportion  to  the  amount  of  this 
interference.   The  presence  of  structural  damping,  however, 
can  partially  alleviate  this  adverse  situation  by  increasing 
flutter  speeds,  provided  the  wing  density  parameter  is  large. 
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Recommendations  for  future  study  include  the  development 
of  the  elementary  theory  for  bending  oscillations,  and  the 
generation  of  numerical  results  for  two-degree-of -freedom 
flutter  using  Program  A.   Of  primary  importance,  however,  is 
the  development  of  a  method  of  characteristics  approach  for 
cascades  with  subsonic  leading-edge  locus.   This  is  the  criti- 
cal problem  today  in  the  design  of  the  supersonic  compressor. 
With  this  problem  properly  formulated,  the  method  of  charac- 
teristics should  be  able  to  give  rapid,  accurate  flutter 
predictions  for  cascades  oscillating  in  this  complicated 
flow  field. 
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APPENDIX  A 
'  COMPUTER  PROGRAM  VARIABLES 

I.   LOGIC  VARIABLES 

The  following  variables  and  their  definitions,  in  alpha- 
betical order  appear  in  the  logic  statements  of  the  programs. 
The  asterisk  is  placed  behind  those  variables  that  appear 
only  in  Program  A. 

FCONST*      Denotes  type  of  oscillation:  TRUE  -  plunge 
mode,  FALSE  -  pitch  mode. 

IATNWL       When  computation  reaches  upper  airfoil,  it 
takes  on  the  value  "1",  otherwise,  "0". 

ICO  Takes  on  the  value  "1"  when  lower  airfoil 

computation  is  complete,  otherwise  "0". 

ICOUNT       Incremented  only  when  MACHLN   is  called  by 

being  set  equal  to  JLINE,  zeroed  when  HIFOIL 
is  called. 

IFIN  Set  equal  to  "1"  when  integration  on  lower 

airfoil  is  complete,  otherwise  "0". 

IHAVEP       The  number  of  the  computation  in  any  row 
starting  with  "0". 

IJUNC        The  grid  point  at  the  intersection  of  the 
initial  right  and  left-running  Mach  lines. 

IJUNT        Equal  to  IJUNC,  needed  to  preserve  the  value 
initially  calculated  in  case  of  multiple 
flow  field  calculations. 

ILINE        The  number  of  grid  points  in  the  row  (right- 
running  Mach  line)  being  computed,  plus  one. 

IMAGRT*      If  X  is  negative,  set  equal  to  "1",  other- 
Wise  "0M. 

IREF         Set  equal  to  "1"  if  lower  airfoil  point  is 
in  Zone  I,  otherwise  "0". 

ISPLIT       Set  equal  to  1/2  of  the  number  of  grid  points 
minus  the  number  of  stagger  points. 
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ISWICH 


ISWTCH 


ITEM 


Set  equal  to  "1"  or  "2"  alternately  starting 
with.  "1"  and  changing  when  upper  airfoil  cal- 
culation goes  into  a  different  zone.   The 
upper  airfoil  zone  marker. 

Switch  variable  used  to  choose  the  row  being 
calculated,  alternately  set  equal  to  "1"  or 
"2".   Used  as  a  second  subscript. 

Loop  variable  for  pitch  and  plunge  mode  cal- 
culation . 


IZOT* 


JCOUNT 


JFIN 


JLINE 


Set  equal  to  "1"  if  flutter  frequency  ratio 
is  imaginary,  otherwise  "0". 

Set  equal  to  "0"  when  in  Zone  I,  incremented 
when  MACHLN  is  called  in  an  interference 
region . 

Set  equal  to  "1"  when  integration  along 
upper  airfoil  is  complete,  otherwise  "0". 

The  number  of  grid  points  in  the  row  being 
computed . 


JSWICH 


The  same  as  ISWICH  but  used  in  LOFOIL.   The 
lower  airfoil  zone  marker. 


JSWTCH 


JTEM 


KCOUNT 


KOUNT 


Switch  variable  used  to  choose  the  preceding 
line  that  was  calculated,  set  equal  to  the 
opposite  value  as  that  of  ISWTCH.   Used  as 
a  second  subscript. 

Loop  variable  for  frequency  range  increment- 
ing. 

Used  in  the  interference  calculation  in 
MACHLN,  set  equal  to  "JCOUNT-1" . 

Set  equal  to  "0"  until  HIFOIL  is  called,  then 
set  equal  to  the  number  of  grid  points  minus 
one . 


LCOUNT 


The  number  starting  with  "0"  of  each  point 
in  a  zone  of  interference  in  LOFOIL. 


MCOUNT 


Same  as  LCOUNT  but  used  in  HIFOIL. 
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II.   QUANTITY  VARIABLES 

The  following  variables  in  the  programs  take  on  the 
values  as  defined  below.   All  have  real  values  unless  other- 
wise stated.   The  asterisk  is  placed  behind  those  variables 
that  appear  only  in  Program  A  and  the  accent  is  placed  behind 
those  variables  that  appear  only  in  Program  B.   The  dimen- 
sioned variables  are  listed  first. 

U22R,  U22I,  V22R,  V22I,  C22R,  C22I(400,3).   The  velo- 
cities, real  and  imaginary,  of  a  point  in 
the  program.   They  are  equivalent  to  the 
respective  Teipel  amplitude  functions, 

U22R  throuSh  C22I' 

U33R,  U33I,  V33R,  V33I,  C33R,  C33I(200,2).   Same  as 

above,  but  used  to  calculate  the  flow  field 
in  Zone  I'  when  the  phase  angle  is  non-zero. 

X,  Y(400,2)  The  x  and  y  coordinates  of  a  point  based  on 
the  chord  length. 

AI  A   as  defined  in  Eq.  (96a). 

BI  B   as  defined  in  Eq .  (96b). 

CAPX'         X  as  defined  in  the  Table  of  Symbols. 

CPU         The  imaginary  part  of  e"1  (c22R+iC22l^ of  the 
preceding  point  calculated  on  the 

upper  airfoil. 

CPR1         The  real  part  of  e"1  (C0 0_+iC0 OT)  of  the 

ing  point  calcu- 
lated on  the  upper  airfoil. 

CPI2         The  imaginary  part  of  e~   (C22  +iC22  )of  the 
point  being  calculated  on  the 

upper  airfoil. 

CPR2         The  real  part  of  e""1  CC£2  tiC22  )  of  the  point 
being  calculated  on  the  upper  airfoil 

DAMP'        The  total  system  damping  as  defined  in  Eq. 
(218b) . 

DELCI        The  value  of  C22   at  x  =  1. 
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DELCPI       The  imaginary  part  of  e~   (DELCR  +iDELCI)  on 
the  upper  airfoil. 

DELCPR       The  real  part  of  e   Q  (JDELCR  +  iDELCI)  on  the 
upper  airfoil. 

DELCR        The  value  of  C  2   at  x  «  1. 

DELL1,  DELL2,  DELM1,  DELM2 .   The  change  in  lift  (real 

and  imaginary)  and  moment,  respectively,  due 
to  extrapolation  to  x  =  1. 

DELTA        The  phase  angle  in  radians. 

DELTAS       The  size  of  the  length  increment  along  a 
Mach  line. 

DISC*        The  discriminant  of  Eq.  (206). 

DSTSTR       The  size  of  the  Ax  increment. 

FACTOR'       Set  equal  to  Ax/k2 . 

FACTOR*      In  the  MAIN,  set  equal  to  Ax/k2  for  the 

pitch  mode,  and  set  equal  to  1/2  Ax/k2  for 

the  plunge  mode.   In  FLUTER,  set  equal  to 

DISC  -  C_. 
K 

FAZE         The  phase  angle  in  degrees. 

FIN  The  chord  length. 

FMANGL       The  Mach  angle  in  radians. 

FNGDPT       The  real  number  value  of  the  number  of  grid 
points  on  the  initial  Mach  line. 

FNGRDN       The  real  number  value  of  the  grid  fineness 
ratio,  equal  to  FNGDPT. 

FSTRMN       The  freestream  Mach  number. 

G1,...,G9    Intermediate  factors  used  in  velocity  calcu- 
lations in  MACHLN,  L0F0IL,  GENFPT,  and  HIFOIL 

GSUBA*       The  structural  damping  coefficient,  g,  as 
defined  in  the  Table  of  Symbols. 

HAFREQ*      The  ratio  of  the  natural  frequency  of  the 
airfoil  in  plunge  to  that  in  pitch. 

HDSTRL       The  increment,  Ax/2. 
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INCRE        The  amount  that  2/k  is  incremented  each  time 
the  lift  and  moment  is  calculated. 

K12R, . . . ,K56I       K12R' *  *  *  ,K56I  aS  defined  in  Ecl  *  <95)- 

LI,  L2       The  real  and  imaginary  lift  coefficients, 
respectively.   They  become  L..  and  L„  or  L- 
and  L,  as  defined  in  Eq . 

(186)and  Eq.  (187)  depending  on  the  type  of 
oscillation . 

L3,  L4*      The  real  and  imaginary  lift  coefficients, 

respectively,  due  to  pitching  oscillations. 

Ml,  M2       The  real  and  imaginary  moment  coefficients, 
respectively.   They  become  M-.  and  M2  or  Mo 
and  Ma  as  defined  in  Eq. 

(188)and  Eq .  (189)  depending  on  the  type  of 
oscillation . 

M3 ,  M4*      The  real  and  imaginary  moment  coefficients, 
respectively,  due  to  pitching  oscillations. 

MUU  The  wing  density  parameter  as  defined  in  the 

Table  of  Symbols. 

NGDPTS       The  grid  fineness  ratio  plus  one. 

NGRDFN       The  grid  fineness  ratio,  it  must  be  an  even 
number  less  than  400.   It  is  an  integer. 

NSTPTS       The  number  of  stagger  grid  points,  an  integer, 
set  equal  to  B/Ax  plus  one. 

OMEGA        The  non-dimensional  flutter  frequency  as 
defined  in  Eq .  (208) . 

OMEGAA  ftfl  as  defined  in  the  Table  of  Symbols. 

OMEGAH*  ft   as  defined  in  the  Table  of  Symbols. 

REDFRQ  The  reduced  frequency. 

RFREQ  2/k. 

R00TX  The  square  root  of  X. 

R00TX1*  The  square  root  of  the  first  root  of  Eq.(206). 

R00TX2*      The  square  root  of  the  second  root  of  Eq. 
(206) . 

ROOTXI*      The  square  root  of  Eq.  (207). 
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RSUBA        The  radius  of  gyration,  rQ,  as  defined  in 
the  Table  of  Symbols. 

RTOSPH  The  ratio  of  specific  heats. 

S  The  square  root  of  M2-l. 

STGANG  The  stagger  angle  in  degrees. 

STGR  The  amount  of  stagger,  B. 

T  M2-l. 

TEST1,  TEST2*   The  difference  between  the  respective 
roots  of  Eq.  (206)  and  that  of  Eq .  (207). 

TEST3, .TEST4*   The  values  of  TEST1  and  TEST2,  respective- 
ly, for  the  previous  frequency  computation. 

The  perpendicular  distance  between  the  two 
airfoils  ,  d  . 

The  chord  length  of  the  top  airfoil,  x-B. 

The  increment,  yAx  cos  a. 

M2 
The  factor,  cos(k   x) . 

M2-l 

The  non-dimensional  flutter  speed  as  defined 

in  Eq.  (209)  . 

UR,  01,  VR,  VI,  CR,  CI   The  values  of  U    , . . . , C    ,  res- 
pectively, for  calculations 

on  the  initial  Mach  line  between  the  point  of 
intersection  of  the  initial  left  and  right- 
running  Mach  lines  and  the  upper  airfoil. 

M2 
V  The  factor,  sin  (k  x) . 

M2-l 
VIPANL       The  imaginary  part  of  the  normal  flow  velo- 
city at  the  airfoil  surface. 

VRPANL       The  real  part  of  the  normal  flow  velocity  at 
the  airfoil  surface. 

W  M2. 

XI*       -  X   as  defined  in  Eq.  (207). 

XLENGTH  The  length  of  the  initial  Mach  line. 

XR1,  XR2*  The  roots  of  Eq.  (206). 
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TNWDST 

TOPCRD 

TRNGLH 

XSUBA*       The  distance  from  the  elastic  axis  of  the 
airfoil  to  its  center  of  gravity. 

XSUBO        The  distance  from  the  leading-edge  of  the 
airfoil  to  its  elastic  axis. 
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APPENDIX  B 


FLOW  DIAGRAM  -  PROGRAM  A 


I.   MAIN  PROGRAM 


Read/Write  DATE 
CALL  INPUT 
IZOT  =  0 


/TToo V 

>\KTEM_=  1,4    / 


^/T0b5        \ 
^v  JTEM  =  1,50  / 


FCONST  =  .FALSE. 


-/TooT         \ 

r  ~^y  ITEM  =  1,2  / 


ICO  =  0 

KOUNT  =  0 

CALL  INITIAL 

1000 

CALL  COMPXY 


IHAVEP.EQ.O 
,AND. 


IATNWL.EQ.O  yT 


IHAVEP.EQ.O 

.AND. 
IATNWL.NE.O 


^IATNWL.EQ.l 

.AND. 
sJHAVEP.EQ.  KOUNT 
"  F 

, 1 s 

(lHAVEP.EQ.JLINE   J- 


J1 


CALL  MACHLN 


CALL   HIFOIL 


CALL  LOFOIL 
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A  ., 


CALL   GENFPT 


IATNWL.EQ.l 

.AND. 
JHAVEP.EQ.KOUNT 


(iHAVEP.NE.JLINE  > 


((X(IHAVEP+1,ISWTCH)+DSTSTR).GT.FIN    J ■»( JLINE- 1 .  LT .  I JUNC ")— ■ 


(iLINE.EQ. NGDPTS > 


ILINE  =   ILINE+1 
JLINE  =   JLINE+1 


ICO  =1 


IATNWL  =  1 
JLINE  =  NGDPTS 
KOUNT  =  NGDPTS- 1 


(iLINE.EQ. NGDPTS > 


IATNWL  =  1 


(iCO.EQ.Q 


I JUNC  =  I JUNC- 1 


Ciswtch.eqTT)- 


ISWTCH  =  1 
JSWTCH  =  2 
IHAVEP  =  0 


ISWTCH  =  2 
JSWTCH  =  1 
IHAVEP  =  0 


(.  NOT.  F  CONST)— 

TXT 


Complete 
L1,L2,M1,M2 


Complete 

L3,L4,M3,M4 


FCONST  =  .TRUE. 


1002 


IHAVEP  =   IHAVEP+1 


ICO.EQ.l 
.AND. 

ijunc.eq.o 
Tf 


ICO.EQ.l 

.AND. 
IHAVEP. GT. I JUNC 
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,.,11 


CALL  FLUTER 


(iMAGRT.NE.l)- 


IZOT  =  1 


RFREQ  =  RFREQ+INCRE 
REDFRQ  =  2 /RFREQ 


1005 


Write  OMEGA,  UF,  JTEM 
XSUBO  =  XSUBO  +0.1 


1100 


STOP 
END 


(TEST3.GT.1E-05) 


, £_ 

< — CTEST4.GT.1E-05") 
T  


f,< 


RFREQ  =  RFREQ- 1 
INCRE  =  1NCRE/5 


-(INCRE.LT.  1E-05) 


JTEM.NE.l 

.OR. 
IZOT.NE.l 


Define  TEST3, 
TEST4 


Define  TESTl, 
TEST  2 


ABS (TESTl ).LE.lE-05 

.OR. 
ABS(TEST2).LE.lE-05 


(TEST1.LT.-1E-05) 


(TEST3.LT. -1E-05 


(TEST2.LT.-1E-05  ) 


(TEST4.LT.-1E-05  >-> 
>j  F  T 


TEST3   =  TESTl 
TEST4   =  TEST2 
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II.      SUBROUTINE   INPUT 


Read /Write  NAMELIST 
Calculate  DSTSTR,  NSTPTS 


(STGR-DTSTSR.  GE.  0  TWy 


(NGDPTS/2.  EQ. NSTPTS- 1  > 


(nGDPTS/2.EQ.  NSTPTS+1  )•; 


*-  NSTPTS  -  NSTPTS+1 


HNSTPTS   =  NSTPTS-1 

l 


NSTPTS  =  NSTPTS+1 

I 


Calculate  STGR,  STGANG,IJUNT, NSTPTS, 

ISPLIT,  S,  DELTA,  OMEGAA,  OMEGAH 

Write  NAMELIST 

RETURN 

END 


III.   SUBROUTINE  INTIAL 


Define  IJUNC,  AI,  BI 


(.'not.fconst  y 


Define  VRPANL,VIPANL 


Define  VRPANL,VIPANL 


Calculate  U22R~>  C22I(1,1)  ,X(1,1) ,  Y(l,l) 


ISPLIT. NE.-l 

.AND. 
DELTA. LE.1E-05 


Calculate  Ll,L2,  M1.M2 


Calculate  Ll,L2,  Ml,M2 


64 
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G3 


1 


Cdelta.le.ie-05  > 


(.NOT.FCONST) 


[Define  VRPANL,VIPANL 


Define  VRPANL,VIPANL 


Calculate  U33R  -•  C33I(1,1) 


(lSPLIT.EQ.-l> 


Calculate  Ll,L2,  M1.M2 


ISWTCH  =  2 
JSWTCH  =  1 
ISWICH  =  1 
JSWICH  =  1 
IFIN  =  0 
JFIN  =  0 
ILINE  =  2 
JLINE  =  1 
IHAVEP  =  0 
IATNWL  =  0 
JCOUOT  =  0 
LCOUNT  =  0 
MCOUNT  =  1 
RETURN 
END 
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IV.   SUBROUTINE  MACHLN 


ICOUNT  =  JLINE 


ICOUNT. LT. I JUNC  ) 


JCOUNT  =  JCOUNT+1 
KCOUNT  =  KCOUNT-1 


(kcount.ne.q) 


Define  UR,UI,  VR,  VI,  CR,  CI 


Calculate  U22R  -  C22I(1,ISWTCH) 


Define  W,T 


I 


(,NOTtFCONST> 


Define  VRPANL,VIPANL 


Calculate  U22R  ~*   C22I(1,ISWTCH) 


XOUNT.NE.ISPLIT 
.OR. 

DELTA. GT.1E-05 


Define  UR,UI,  VR,  VI,  CR,  CI 


J 


Define  VRPANL,VIPANL 


)ELTA.LE.lE-05 

.OR. 
JLINE.  GT.  IS PLIT 


Define  U22R  ->  C22I(1,3) 


■(.NOT.FCONST  ) 
I  F 


Define  VRPANL,  VIPANL 


Define  VRPANL,VIPANL 


I 


(ICOUNT. NE. IS PLIT  )- 


Define  U22R  ->  C22I(1,3) 


Calculate  U33R  -  C33I(1, ISWTCH) 


RETURN 
END 
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V.   SUBROUTINE  HIFOIL 


Define  TOPCRD 


(JLINE.  EQ.  NGDPTS-1*)- 


IATNWL  =  0 


(ICOUNT.EQ.O  > 


(JCOUNT.NE.O> 


Define  K12R,K12I,J 


Define  K12R,K12I,J 


([COUNT.  NeTQ^)- 


Define  I,  J,  K12R,K12I 


(.NOT.FCONST  )- 


Define  V22R,V22I(1,ISWTCH) 


1  Define  V22R,V22I(1,ISWTCH) 


Calculate  U22R,  U22I,  C22R,  C22I(1,ISWTCH) , 
CPR1,  CPU,  CPR2,  CPI2 


(jFIN.EQ.Q 


.COUNT. EQ.O 

.AND. 
?OPCRD+DSTSTR. GT . 1 . 0 


'LCOUNT.EQ.O 
.OR. 
:OPCRD+DSTSTR. GT . 1 . 0 


Calculate  Ll,L2,Ml,M2 


Calculate  Ll,L2,Ml,M2 


LCOUNT   =  LCOUNT4-1 


130 
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( ISWICH. NE.1> 


(  LCOUNT.  EQ.IJUNT). 


ISWICH  =  2 


(  JCOUNT.EQ.O): 


->-«- 


[m 


(  LCOUNT . NE . JCOUNTp- 
t\    F 


Calculate  Ll,  L2,  Ml,M2 
LCOUNT  =  0 


(  ISWICH.  EQ.L) 


ISWICH  =  1 


(TOPCRD+DSTSTR.LE.1.0  )- 


Calculate  Ll,  L2,  Ml,  M2 
JFIN  =  1 


ICOUNT  =  0 

RETURN 

END 
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VI.       SUBROUTINE   GENFPT 


JDefine   I 


(iatnwl.eqTcT)- 


Define  K12R,   K12I,    K34R,    K34I  [  [Define  K12R,    K12I,   K34R,   K34I 


Calculate  U22R  -»  C22I(I,ISWTCH) 


1COUNT.NE.ISPLIT 
.OR. 

DELTA. GT. IE- 05 


Define, U22R  ■*  C22I(I,3) 


DELTA. LE.1E-05 
.OR. 

JLINE.GT.ISPLIT 


Define  K12R,   K12I,    K34R,   K34I,    K56R,   K56I 
Calculate  U33R  -  C33I(I,ISWTCH) 


(JLINE.NE.ISPLIT> 


Define  U22R  -»     C22I(I,3) 


RETURN 
END 
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vii.  subroutine  lofoil 
(.not.fconst)— 


Define  V22R,V22I(ILINE,ISWTCH)|    | Define  V22R,V22I (ILINE, ISWTCH) 


Define  K12R,  K12I.  K34R.  K34I 


(  IATNWL.EQ.O> 


Define  K56R,  K56I 


Define  K56R,  K56I 


Calculate  U22R,U22I, C22R, C22I (ILINE, ISWTCH) 


c 


ICOUNT.NE.ISPLIT 
.OR. 

DELTA. GT.1E-05     /  T 


Define  U22R  ->  C22I(ILINE,3) 


DELTA. LE.lE-05 
.OR. 
T  V  JLINE.GT.ISPLIT 


(.NOT.FCONST  ) 


Define  V33R,V33I (ILINE, ISWTCH)     Define  V33R,V33I (ILINE, ISWTCH) 


Calculate  U33R,U33I,C33R,C33I (ILINE, ISWTCH) 


(JLINE.NE.ISPLIT    > 


Define  U22R  -»  C22I(ILINE,3) 


88 


■>-  94 
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Define   CPRl,    CPU,    CPR2,CPI2 


(IFIN.EQ.T) 


Calculate  Ll,L2,   Ml,   M2 


JLINE.NE.ISPLIT 

.AND. 
X(K,L)+DSTSTR.LE.1.0 


Calculate  Ll,L2,  Ml,  M2 


IREF  =  0 


ClFIN.EQ.l  > 


'MCOUNT.EQ.O 

.AND. 
X(ILINE, ISWTCH)+DSTSTR.GT. 1 .0 


140 


MCOUNT.EQ.O 

.OR. 
,X(ILINE,ISWTCH)+DSTSTR.GT.1.0  ^/T 
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Calculate  Ll,L2,  Ml,  M2 


Calculate  L1,L2,  Ml,  M2 


(iREF.EQ.T") H  130 


<IREF.EQ.l) 


MCOUNT  =  MCOUNT+1 


(X(ILINE,ISWTCH)+DSTSTR.GT.T70 


(jSWICH.NE.l> 


110 


" 1 v 

(mcount.ne.  isplit+Q- 


(mcount.eq.ijunc    y 


Calculate  Ll,L2,   Ml,   M2 


I 


1  120k 


-(iREF.EQ.l) 

irpr. 
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^ 

i  /. 

Va(.ILIin£j,  j.swnJH.;-t-uaxbi"JK..iJJE..  i  .  u 

)T 

14-u 

1 

F 

\ 

f 

Calculate  DELLl,DELL2,DELMl,DELM2 

Calculate  Ll,L2,   Ml,   M2 

' 

' 

(DELT 

A    CT    lE-n5^                 ■*. 

«■» 

F                      T 

' 

1 

Calculate  Ll,L2,   Ml,   M2 
IFIN  =  1 

ClREF 

.  LQ.UJ                               ^ 

F                      T 

_J 

.s 

\ 

■ 

Cmcou 

TTT    PT    TnPT  TT  1  1  J  . ..     — 

RETURN 

lUi  LtJ-.  i.  i_>  1  J_i  1.  J.  1  L/ 

F                                  T 

ENI 

) 
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VIII.   SUBROUTINE  COMPXY 


(IHAVEP.EQ.0> 


Define  I 
Calculate  X,Y(I,ISWTCH) 

zi — z= 


Define  I 

Calculate  X,Y(I,ISWTCH) 


RETURN 
END 


IATNWL.EQ.l 

.AND. 
JLINE.NE.NGDPTS-1 


Define  I 

Calculate  X,Y(I,ISWTCH) 


IX.   SUBROUTINE  FLUTER 


IMAGRT  =  0 
Calculate  CR,CI 


X61 


OMEGAH.LT.1E-06 
.OR. 

OMEGAA.LT.1E-06 


[Calculate  DISC 


(disc.le.-ie-Q6  J- 


|  Calculate  XR1,XR2~ 


1  Calculate  XI 


XR1.LE.-1E-06 

.OR. 
XR2.LE.-1E-06 

.OR. 
XI.LE.-1E-06 


Calculate  OMEGA, UF 

Write  RFREQ,R00TX1,R00TX2,R00TXI 


RETURN 
END 


i 


Calculate  XRl,XR2 


IMAGRT  =   1 
Write  RFREQ 
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